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evaluation of pulmonary nodules from Computed Tomography (CT) scans: detection, in which the locations of possible nodules are 
identified, and characterization, in which a nodule is represented by measured features that may be used to evaluate the probability 
that the nodule is cancer. Currently, the most useful prediction feature is growth rate, which requires the comparison of size estimates 
from two CT scans recorded at different times. The present invention includes methods for detection and feature extraction for size 
characterization. The invention focuses the analysis of small pulmonary nodules that are less than 1 centimeter in size, but is also 
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SYSTEM. METHOD AND APPARATUS FOR SMALL PULMONARY 
NODULE COMPUTER AIDED DIAGNOSIS FROM 
COMPUTED TOMOGRAPHY SCANS 

This application claims the benefit of U.S. Provisional Application 
No.60/322,038, filed September 14, 2001, which is incorporated herein by 
reference. 

BACKGROUND OF THE INVENTION 

5 The present invention relates to the art of diagnostic imaging of small 

pulmonary nodules. In particular, the present invention is related to analyzing and 
manipulating computed tomography scans to: segment the lungs, measure lung 
volume, locate and determine the size of the nodules without explicit segmentation, 
register the nodules using a rigid-body transformation, and removing the pleural 
1 0 sxuf ace firom juxtapleural nodules in thresholded images. 

Lung cancer is the leading cause of cancer deaths among the population in 
the United States. Each year there are about 170,000 newly diagnosed cases of 
lung cancer and over 150,000 deaths. More people die of lung cancer than of 

15 colon, breast, and prostate cancers combined. Despite the research and 

improvements in medical treatments related to surgery, radiation therapy, and 
chemotherapy, currently the overall survival rate of all lung cancer patients is only 
about 14 percent. Unfortunately the survival rate has remained essentially the 
same over the past three decades. The high mortality rate of lung cancer is caused 

20 by the fact that more than 80% lung cancer is diagnosed after it has metastasized. 
Patients with early detection of lung cancer followed by proper treatment with 
surgery and/or combined with radiation and chemotherapy can improve their five- 
year survival rate from 13 percent to about 41 percent. Given that earlier-stage 
intervention leads to substantially higher rates of survival, it is therefore a major 

25 public health directive to reduce the mortality of lung cancer through detection and 
intervention of the cancer at earlier and more curable stages. 
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The development of the computed tomography (CT) technology and post- 
processing algorithms has provided radiologists with a useful tool for diagnosing 
lung cancers at early stages. However, current CT systems have their inherent 
shortcomings in that the amount of chest CT images (data) that is generated from a 
5 single CT examination, which can range from 30 to over 300 slices depending on 
image resolution along the scan axial direction, becomes a huge hurdle for the 
radiologists to interpret. Accordingly, there is a constant need for the improvement 
and development of diagnostic tools for enabling a radiologist to review and 
interpret the vast amount of information that is obtained through a CT examination. 

10 

International Publication No. WO 01/78005 A2 discloses a system and 
method for three dimensional image rendering and analysis, and is incorporated 
herein by reference. The system performs a variety of tasks that aid a radiologist 
in interpreting the results of a CT examination. 

15 

One task that radiologists focus on is segmenting the lung region from the 
image of a single slice obtained from the CT examination, hi the prior art, some 
have suggested using a linear discriminant function and morphological filtering to 
automatically segment the lungs (S. Hu, E. A. Hoffman, and J. M. Reinhardt, 

20 "Automatic Lung Segmentation for Accurate Quantitation of Volumetric X-Ray 
CT Images," IEEE Transactions on Medical Imaging, Vol 20, No 6, June 2001, 
which is incorporated herein by reference) while others have used mean and 
median filtering to remove the streaking artifacts due to excessive x-ray quantum 
noise (J. Hsieh, "Generalized Adaptive Median Filters and their Application in 

25 CT," SPIE, Vol 2299, 1994; J. Hsieh, "Adaptive Trimmed Mean Filter for CT 
Imaging," SPIE, Vol 2299, 1994 which is incorporated herein by reference). 

Radiologists also study the location and size of the pulmonary nodules in 
the CT scan. It is preferred if the radiologist could perform this analysis without 
30 the use of explicit segmentation. In some prior work, the location of a nodule was 
determined by finding the center of mass of the nodule through an iterative 
correlation-based procedure (A.P. Reeves, W.J. Kostis, D.F. Yankelevitz, CI. 
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Henschke, "Analysis of Small Pulmonary Nodules without Explicit Segmentation 
of CT images," Radiological Society of North America - 2000 Scientific Program, 
vol. 217, pgs. 243-4, November 2000 which is incorporated herein by reference). 
The method works for isolated pulmonary nodules, but fails on nodules attached to 
5 the pleural surface. 

Radiologists also estimate a measurement of doubling time of a nodule by 
registering two separate images of the nodule taken at two different times (time-1 
and time-2). This analysis requires that the time-1 and time-2 nodules be 

10 registered correctly so that the growth can be properly measured. Other objects 
such as vessels and bronchial tubes must also be registered together. This results 
in their absence in the difference unage and little effect on the growth 
measurement. Previously, two nodules were registered by finding the centers of 
mass of the nodules and translating the image accordingly (A. P. Reeves, W. J. 

1 5 Kostis, D. F. Yankelevitz, C. I. Henschke, "Analysis of Small Pulmonary Nodules 
without Explicit Segmentation of CT images " Radiological Society of North 
America - 2000 Scientific Program, vol. 217, pgs. 243-4, November 2000 which 
is incorporated herein by reference). However, this analysis did not guarantee that 
the two nodules would be correctly orientated, and that the other objects in the 

20 image would registered because these objects might be rotated about the nodule. 
Some have registered the nodules by performing a maximization search of the 
mutual information metric over the rigid-body transformation parameters (F. Maes, 
A. CoUignon, D. Vandermeulen, G. Marchal and P. Suetens, "Multimodality image 
registration by maximization of mutual information," IEEE Transactions on 

25 Medical Bnaging, vol 16, no. 2, pgs. 187-198, April 1997; Takagi, N.; Kawata, 
Y.; Nikvi, N.; Morit, K.; Ohmatsu, H.; Kakinuma, R.; Eguchi, K.; Kusumoto, M.; 
Kaneko, M.; Moriyama, N. "Computerized characterization of contrast 
enhancement patterns for classifying pulmonary nodules" Image Processing, 2000, 
Proceedings. 2000 International Conference on, vol, l,pgs. 188-191,2000 

30 which are incorporated herein by reference). 



3 



wo 03/024184 



PCT/US02/29366 



Radiologists also need to remove the pleural surface ftom juxtapleural 
nodules in CT images. In some prior work, three-dimensional morphological 
filtering and mathematical moments were used to segment a juxtapleural nodule 
from pleural surface in a binary image (A. P. Reeves, W. J. Kostis, "Computer- 
5 Aided Diagnosis of Small Pulmonary Nodules," Seminars in Ultrasound, CT, and 
Mi?/, vol. 21, no. 2,pgs. 116-128, April 2000 which is incorporated herein by 
reference). 

SUMMARY OF THE INVENTION 

10 The present invention is directed to diagnostic imaging of small pulmonary 

nodules. There are two main stages in the evaluation of pulmonary nodules from 
Computed Tomography (CT) scans: detection, in which the locations of possible 
nodules are identified, and characterization, in which a nodule is represented by 
measured features that may be used to evaluate the probability that the nodule is 

1 5 cancer. Currently, the most useful prediction feature is growth rate, which requires 
the comparison of size estimates from two CT scans recorded at different times. 
The present invention includes methods for detection and feature extraction for 
size characterization. The invention focuses the analysis of small pulmonary 
nodules that are less than 1 centimeter in size, but is also suitable for larger nodules 

20 as well. 

For the purpose of Computer Aided Diagnosis (CAD), pulmonary nodules 
are dichotomized into attached nodules and isolated nodules based on their location 
with respect to other solid lung structures. Attached nodules are adjacent to some 
25 larger solid structure, such as the pleural surface. Isolated nodules consist of both 
well-circumscribed nodules and nodules that are larger than all adjacent structures, 
such as blood vessels or bronchi. Nodules themselves may be solid, non-solid or 
part-solid. The analysis of a CT scan for the existence and study of pulmonary 
nodules generally entails the following: 

30 
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1 . Detection 

(a) Identify the lung regions and main bronchi from thoracic CT images 

5 (b) Separate the lungs into two major regions: (1) the lung parenchyma 

and (2) the lung surface region, including the pleural surface and major airways. 

(c) Identify possible locations of isolated nodules in the lung 
parenchyma region and identify possible locations of attached nodules in the in the 
1 0 lung surface regions. 

2. Characterization 

(a) Starting with a single location point within a possible nodule, 

15 identify the nodule region in the CT images. This entails locating the geometric 
center of the nodule and approximating its size. 

(b) Given the location and approximate size of a nodule, compute 
characteristic features of the nodule, including robust size estimates. 

20 

In connection with the overall methodology of analyzing CT scans, the 
present invention includes sub-methods for the following: 

1 . The segmentation of whole lung CT scans into lung parenchyma 
25 and lung surface regions 

The first preprocessing stage in CT lung image analysis is to obtain the 
regions of interest from the whole lung scans. The lung region consists of all 
tissue found within the pleural surface, including lung parench3mia, vessels, and 
30 possibly nodules. Features of this approach are the partitioning of the lung into 

three major regions and the tailoring of the segmentation algorithm for each region. 
In addition a distinction is made between the central lung parenchyma and the 

5 
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region of the of the lung parenchyma near to the lung walls. A properly segmented 
lung region greatly reduces the search space of an entire CT scan. 

2. The characterization of lung air volume and inspiration both from 
5 the entire parenchyma region and from single axial CT images 

The total volume of the lungs is estimated from the whole lung scans. In 
addition the change in volume between different scans due to changes in 
inspiration is also addressed. 

10 

3. The automatic location and size characterization of nodules 

An algorithm finds the center and approximate size of pulmonary nodules 
in CT images without the use of explicit segmentation. The algorithm has a 
1 5 weighting function that locates the center of mass of the nodule. A second 
weighting ftmction, centered on the location of the first weighting ftinction, 
estimates the size of the nodule. The process is repeated with weighting ftmctions 
of increasing size until the size of the nodule region is reliably determined. The 
algorithm works for both isolated nodules and nodules attached to the walls. 

20 

The algorithm is designed for use on images that are resampled from high- 
resolution CT scans (1mm slice thickness) into an isotropic voxel space. Such a 
system could be useful in helping radiologists analyze pulmonary nodules. In the 
CT image browser, the radiologist could click on a nodule and, using the 
25 algorithm, the computer could automatically locate the center and size of the 

nodule. Using this information, the computer then clips out a region of interest for 
the nodule and performs some nodule characterization analysis and perhaps some 
3D visualization. 

30 4. The registration of nodules from two scans 
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An algorithm registers two different scans of a nodule region. A three- 
dimensional rigid-body transformation of one scan is made to optimally match the 
location and orientation (6-degrees of freedom) of the second scan. Powell's 
method is the preferred search strategy. 

5 

The analysis of pulmonary nodules without explicit segmentation estimates 
a measurement of doubling time by applying a weighted function on the difference 
image of the two nodules. The analysis requires that the time-1 and time-2 nodules 
be registered correctly so that the growth can be properly measured. Other objects 
1 0 such as vessels and bronchial tubes must also be registered together, resulting in 
their absence in the difference image and little effect on the growth measurement. 

5. The segmentation of nodules attached to the pleural surface. 

1 5 An algorithm segments nodules that are attached to the pleural surface from 

the pleural surface. Starting from a location within a nodule, the direction of the 
pleural surface is determined. A cutting plane is then iteratively moved towards 
the surface until the volume of the region behind the wall suddenly increases. The 
increase in volume indicates that the pleural surface has been reached. The 

20 orientation and location of the plane is then preferably optimized by a hill climbing 
procedure. 



Tlie algorithm removes the pleural surface (and any other extraneous 
objects) from a high-resolution image containing a juxtapleural nodule. The 

25 algorithm is not only required to perform a good segmentation, but expected to 
perform the same segmentation when given different scans of the same nodule. 
This property is necessary because the resulting segmented nodule will be used for 
the characterization and analysis of the nodule. If the segmentation is not 
performed consistently between images scanned at different times, then the 

30 doubling time estimation will not be accurate. 



7 
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A preferred embodiment of the present invention is a method and apparatus 
for generating a lung mask for segmenting a lung image from voxel data 
containing the lung image, noise, solid components and surrounding background 
information. The apparatus of the invention is a masking unit configured with the 
5 teachings of the method of the invention. The invention can also be practiced on a 
machine readable medium. The method includes initially applying a median filter 
and a mean filter to the voxel data to reduce noise. The noise reduced voxel data is 
then thresholded to identify solid components and other structures. The 
surrounding background in the noise reduced voxel data is identified and deleted. 
10 The connected components of the voxel data are labeled, and the largest connected 
components are determined to select a lung region having a geometric form. The 
voxel data is morphological filtered to refine the geometric form of the lung 
region. 

15 In a preferred embodiment, the voxel data is associated with a plurality of 

slices through a patient's lung with the slices beginning at about the patient's 
shoulders and the application of the median filter and the mean filter is performed 
for the first 25 percent of the plurality of slices only to substantially reduce the 
computation time. Preferably the median filter has a size of 4x4, and the mean 

20 filter has a size of 1x3. Preferably the voxel data in step is thresholded at a gray 

level of about 500. Preferably the largest connected components are selected to be 
associated with more than about 1 percent of the voxel data. Preferably the voxel 
data is associated with a plurality of slices for the morphological filtering with the 
slices being divided into a first end region, a middle region, and a second end 

25 region. The morphological filtering is preferably performed with 2D circular filter 
having a first diameter in the first end region and the second end region while 
being performed with 2D circular filter having a second diameter which is about 
twice the first diameter in the middle region. 

30 The present invention also includes a method and apparatus for measuring 

lung volume fiom a segmented lung image. The lung image is obtained from a 
scan which includes a plurality of slices of voxel data having a gray level value and 

8 
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a volume associated therewith. The apparatus of the invention is a lung volume 
measuring unit configured with the teachings of the method of the invention. The 
invention can also be practiced on a machine readable medium. The method 
includes initially generating a matrix of entries from the segmented lung image. 
5 The matrix includes a plurality of columns and a plurality of rows. Each of the 
plurality of columns represent the gray level value and each of the plurality of rows 
represent one of the plurality of slices in the scan with the entries corresponding to 
the niunber of times that the gray level occurs in the corresponding slice. The 
number of voxels in the segmented lung image is next determined from the matrix 
1 0 entries. The number of voxels in the segmented lung image is multiplied by the 
volume of each voxel. 

The present invention also includes a method and apparatus for measuring 
volume of tissue in a segmented lung image. The lung image is obtained from a 

15 scan which includes a plurality of slices of voxel data having a gray level value and 
a volume associated therewith. The apparatus of the invention is a lung tissue 
measuring unit configured with the teachings of the method of the invention. The 
invention can also be practiced on a machine readable medium. The method 
includes initially generating a matrix of entries from the segmented lung image. 

20 The matrix includes a plurality of colunms and a plurality of rows. Each of the 

plurality of columns represent the gray level value and each of the plurality of rows 
represent one of the plurality of slices in the scan with the entries corresponding to 
the number of times that the gray level occurs in the corresponding slice. The sum 
of tissue of voxels in the segmented lung image is next determined from the matrix 

25 entries. The sum of tissue of voxels in the segmented lung image is multiplied by 
the volume of each voxel. The sum of tissue of voxels is preferably calculated by 
summing the product of each matrix entry multiplied by the corresponding gray 
level value divided by a gray level value assigned for tissue. 

30 The present invention also includes a method and apparatus for finding the 

location, P\ and size, r\ of a pulmonary nodule in a high-resolution computed 
tomography (CT) image. The apparatus of the invention is a nodule finding unit 

9 
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configured with the teachings of the method of the invention. The invention can 
also be practiced on a machine readable medium. In the method, a set of initial 
processing parameters including an initial location, Pu an initial size, ru and target 
value, T, is initially selected. An initial new location i^' of the nodule is computed 
5 with a locator template function. After incrementally increasing size, r/, a new 

location /J' of the nodule is computed with the locator template ftinction, and a size 
metric is computed with a sizing template function. The size, r/, is incrementally 
increased while determining a new location i^/ of the nodule and computing the 

size metric until the size metric is less than the target value. The location, P', and 
10 the size, r', is then retumed j&om the previous iteration of increasing size, 

A preferred method for finding the location, P, and size, r, of a pulmonary 
nodule in a high-resolution computed tomography image initially includes 
windowing the image to ignore bone structures, and selecting a locator template 

1 5 function and a sizing template function. A set of initial processing parameters 
including: an initial location, P; size, r; and termination criteria is selected. A 
search is performed to determine a maximum response of the locator template 
function. A response of the sizing template function is determined and compared 
to the termination criteria. If the termination criteria has not been satisfied, the 

20 size, r, is incrementally increased, the response functions are determined and 

compared to the termination criteria. Once the termination criteria are satisfied the 
location, P, and size, r, of the nodule are outputted. 

In the preferred method for finding the location, P, and size, r, of a 
25 pulmonary nodule in a high-resolution computed tomography image, preferably 
the image at intensities over about 1000 are clipped to window the image. 
Preferably the locator template function is either; a Gaussian template function, a 
Laplacian of the Gaussian template function, or a difference of Gaussians template 
function. The template function preferably has at least four parameters 
30 corresponding to the x-location, y-location, z-location, and radius. The initial 
location, P, can generally be either calculated from the image or specified by a 



10 
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user. Preferred methods for searching include a hill climbing method and Powell's 
method. 

The present invention also includes a method and apparatus for registering 
5 3-d images of a pulmonary nodule from a high-resolution computed tomography 
(CT) scans. The images include a first image (imi) obtained at time-1 and a second 
image (im2) obtained at time-2, and are in a floating point pixel-format associated 
with a 6-dimensional parameter space. The apparatus of the invention is a 
registering unit configured with the teachings of the method of the invention. The 

10 invention can also be practiced on a machine readable medium. The method 

includes calculating initial rigid-body transformation parameters for a rigid-body 
transformation on the first image (zmi). The optimum rigid-body transformation 
parameters are determined by calculating a registration metric between the second 
image (im2) and the rigid-body transformation on the first image (imi). A 

15 registered image is generated from the optimum rigid-body transformation 
parameters. 

In a preferred method for registering 3-d images of a pulmonary nodule 
from a high-resolution computed tomography (CT) scans, the calculation of the 
20 initial rigid-body transformation parameters is preferably preceded by masking one 
of the images by setting pixels to a background value. Preferably the background 
value is about -1000. The registration metric is generally either minimized or 
maximized. In one preferred embodiment, the registration metric is preferably 
calculated by: 

25 transforming the first image (/mi) with the initial rigid-body transformation 

parameters to obtain a transformed first image (imu); 

calculating the registration metric as a correlation (C.) between the 
transformed first image (imu) and the second image (irriz); and 

searching for the maximum correlation (C.) in the 6-dimensional parameter 

30 space. 

In another preferred embodiment, the registration metric is preferably calculated 
by: 

11 
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10 



transforming the first image {im\) with the initial rigid-body transformation 
parameters to obtain a transformed first image (imuy, 

calculating the registration metric as a mean-sqnared-difference (MSD) 
between the transformed first image {im\i) and the second image (/mi); and 

searching for the minimmn mean-squared-difference (MSD) in the 6- 
dimensional parameter space. 

The transforming of the first image {im\) to obtain the transformed first image 
(zmu) is preferably a mapping of a point v in 3-d space to a point V in transformed 
space defined by : 



wherein R^^ Ry^ and i?^ are rotation matrices defined as: 



15 



20 



R. = 



R. = 



10 0 " 

0 cos(r,) -sin(r,) 
0 sin(r^) cos(rJ 

cos(r^) 0 sin(r^) 

0 1 0 
-sin(r_^) 0 cos(r^) 



cos(rJ 
sin(;'L ) 
0 



-sin(;;) 
0 



The initial rigid-body transformation parameters preferably include six parameters 
(tx,ty,tz,rx,ry^rz) respectively defined as translation in x, translation in y, translation 
in z, rotation about the x-axis, rotation about the y-axis, and rotation about the z- 
axis. Preferably the initial rotation parameters (rx,ry,rz) are all set to zero, and the 
initial translation parameters (pc,ty,tz,) are set so that the nodule in the first image 
(imi) overlaps the nodule in the second image (imz) during the initial calculation of 
the registration metric. The initial translation parameters (pc,ty,tz,) can be set to a 
difference between the center of the first image (imi) and the center of the second 
image (imz). Preferably the initial translation parameters (tK,ty,tz) are set to a 

12 
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difference between the center of mass of the first image (imi) and the center of 
mass of the second image (rma). The searching is can be conducted by calculating 
the correlation (C) or mean-squared-difference (MSD) for every possible set of 
rigid-body transformation parameters. Preferably the searching is conducted by 
5 either a Hill-Climbing search method or by Powell's method. 

The present invention also includes a method and apparatus for removing 
extraneous matter fi:om an image having a juxtapleural nodule. The apparatus of 
the invention is a processing unit configured with tihe teachings of the method of 

10 the invention. The invention can also be practiced on a machine readable medium. 
The method includes providing an initial location P'. A spherical volume that fits 
inside the image and is centered at the initial location P' is calculated. A center of 
mass, COM, of the spherical volume is calculated. An initial direction d' directed 
towards the extraneous matter is determined in accordance with the following 

15 equation: 

COM-F 
||COM-P|^ 

A current location P/ is initiaUzed to be equal to the initial location P'. A current 
direction di is initialized to be equal to the initial direction d\ A maximum ratio 
7max> step size s, prior mass massi.u and prior change in mass A/.i are initialized. 

20 The current location Pi is moved by the step size s in the current direction du An 
equation defining a plane A is determined so that the plane A is normal to the 
current direction di and the plane A passes through the current location P/. A 
current mass, massi, is calculated of the nodule on a side of the plane A opposing 
that of the extraneous matter. A current change in mass A/ is calculated by 

25 subtracting the prior mass massi^i from the current mass massi. A current ratio y is 
calculated in accordance with the following equation: 



The prior mass massi.x is set equal to the current mass massi. The prior change in 
mass Af-i is set equal to the current change in mass A/. The current ratio 7 is 
30 compared to the maximum ratio 7max. The current direction di is modified to 

13 
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minimize the current mass massu and the above steps are repeated starting with 
moving the current location P/ by the step size s in the current direction di while the 
current ratio 7 is one of less than and equal to the maximum ratio 7max. The area of 
the nodule partitioned by the plane A is output in response to the current ratio 7 
5 being greater than the maximum ratio 7max. 

In the preferred method for removing extraneous matter from an image 
having a juxtapleural nodule, the extraneous matter can include a pleural surface. 
Preferably the maximum ratio 7max is initialized to 0.5 and initializing the step size 
10 s is initialized to 1 .5. Preferably the method also includes following the steps after 
the current ratio 7 is determined to be greater than the maximum ratio 7max: 
defining the current location Pi as being visited; 

determining on which side of the plane A the current location P, is located; 
terminating in reponse to the current location Pi not being located on a side 
15 of the plane opposing that of the extraneous matter; 

defining the current location P/ as being part of a region of interest in 
response to the current location Pi being located on the side of the plane opposing 
that of the extraneous matter; and 

performing these additional steps recursively using a location 
20 corresponding to at least one of six (6) one-pixel moves from the current location 
P/. 

Preferably the step where the current direction di is modified to minimize the 
current mass massi, and the step where the area of the nodule partitioned by the 
plane A is output in response to the current ratio 7 being greater than the maximum 
25 ratio 7niax include the steps of: 

calculating recursively the current mass massi of the nodule on a side of the 
plane A opposing that of the extraneous matter using at least one of six (6) 
directions and a step size Si from the current location P,; and 

defining the current direction di equal to the direction yielding the largest 
30 decrease in the current mass masSi. 

Preferably the initial location P' is located near a center of the nodule. 
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For a better understanding of the present invention, reference is made to the 
following description to be taken in conjunction with the accompanying drawings 
and its scope will be pointed out in the appended claims. 

5 BRIEF DESCRIPTION OF THE DRAWINGS 

Preferred embodiments of the invention have been chosen for purposes of 
illustration and description and are shown in the accompanying drawings, wherein: 

10 Figure 1 illustrates at the top six (6) images that are selected slices from a 

single computer tomography scan along with the corresponding extracted lung 
regions from the selected slices at the bottom; 

Figure 2 illustrates a histogram of image intensity of the lung region that 
1 5 includes parenchyma, vessels, and nodules and other solid structures such as bone 
and organs; 

Figure 3 is a table illustrating the densities of structures within the lungs; 

20 Figure 4 illustrates an air pocket that is not part of the lungs but has similar 

characteristics to the modeling of the lungs; 

Figure 5 is a lung mask generation algorithm; 

25 Figure 6 is a flow chart of the lung segmentation algorithm with sample 

images incorporated therein; 

Figure 7 illustrates on the left hand side the horizontal streaking artifacts 
due to beam hardening that can be reduced as shown on the right hand side with 
30 median and mean filtering; 
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Figure 8 illustrates an original computer tomography scan on the left hand 
side with the image shown after thresholding in the center and the resulting image 
after the background has been identified, removed and inverted shown on the right 
hand side; 

5 

Figure 9 illustrates on the left hand side a lung after thresholding without 
median and mean filtering and on the right hand side the image after median and 
mean filtering; 

10 Figure 10 illustrates two lungs, gas, and some noise, and that the retention 

of only the largest components would eliminate the unwanted regions; 

Figure 1 1 illustrates at the left hand side a segmented lung after 
thresholding but before morphological filtering, and at the right hand side the lung 
1 5 after morphological closing; 

Figure 12 illustrates the lungs being divided medially into three regions 
where large vessels are found in region 2 which require a large structuring element, 
and that only small vessels are found in regions 1 and 3; 

20 

Figure 13 illustrates a table illustrating the timed execution of 
morphological filtering using both fixed (28 pixel diameter) and varying (28 pixel 
and 14 pixel diameter) structuring elements; 

25 Figure 14 illustrates a matrix containing the histograms of the segmented 

lung slices; 

Figure 15 illustrates an algorithm for finding the location and size of a 

nodule; 

30 

Figure 16 illustrates a one dimensional model of a) isolated nodule next to 
blood vessels, and b) a nodule attached to the pleural wall; 
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Figure 17 illustrates locator template functions; 

Figure 18 illustrates one dimensional examples of the maximizing locations 
5 of template functions for (a-c) a Gaussian template on an isolated nodule, (d-f) 
LOG template on a nodule on the wall, and (g-i) LOG template on an isolated 
nodule; 

Figure 19 illustrates two dimensional examples of the maximizing locations 
10 for two dimensional LOG templates of various radii on isolated and pleural 
nodules where the iimer and outer circles represent the positive and negative 
regions, respectively, of the LOG template; 

Figure 20 illustrates two dimensional examples of a sizing template 
1 5 function of various radii which exhibit dramatic changes in the distribution of 
values inside the circle when the sizing template function becomes too large; 

Figure 21 illustrates a detailed algorithm for finding the location and size of 
a nodule; 

20 

Figure 22 illustrates a response of the mean filter function for various sized 
spheres (solid intensity of 1000) versus the radius of the template function; 

Figure 23 illustrates a response of a Gaussian template function for various 
25 sized spheres (solid intensity of 1000) versus the radius of the template function; 

Figure 24 illustrates a graph of the Gaussian and the Laplacian of the 
Gaussian; 

30 Figure 25 illustrates a hill climbing search algorithm for the nodule finding 

algorithm; 
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Figure 26 illustrates locator template functions; 

Figure 27 illustrates an algorithm for registering a first image and a second 
image using a rigid-body transformation; 

5 

Figure 28 illustrates a hill climbing search algorithm for determining the 
best rigid-body transformation parameters for the registering algorithm; 

Figure 29 illustrates an algorithm for removing the pleural surface from 
1 0 juxtapleural nodules; 

Figure 30 illustrates an example of the pleural wall removal algorithm: The 
cut nodule is shown as dark gray white the rest of the nodule and the pleural 
surface are shown in light gray. (A) Given the initial location and direction, (B) 

15 the plane is moved in the direction. (C) The plane eventually intersects the pleural 
wall causing an increase in the size change of the cut nodule, and (D) the plane is 
reorientated to minimize the size of the cut nodule. At the next iteration, (E) the 
plane intersects the pleural wall, but this time (F) reorientating the plane to 
minimize the cut nodule still results in a large increase in size change. The 

20 algorithm terminates, returning the cut nodule from the previous iteration, (D); 

Figure 31 illustrates an algorithm for recursively finding the cut nodule 
region; and 

25 Figure 32 illustrates a hill climbing search algorithm for the algorithm for 

removing the pleural surface. 

DF.TATT.TCD DESCRIPTION OF THE INVENTION 

30 A system in accordance with the present invention may include a scaruier, 

processor, memory, display device, input devices, such as a mouse and keyboard, 
and a bus coimecting the various components together. The system may be 
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coupled to a communication mediimi, such as a modem coimected to a phone 
line, wireless network, or the Internet, 

The present invention is preferably implemented using a general purpose 
5 digital computer, microprocessor, microcontroller, or digital signal processor 
programmed in accordance with the teachings of the present specification, as will 
be apparent to those skilled in the computer art. Appropriate software coding may 
be readily be prepared by skilled programmers based on the teachings of the 
present disclosure, as will be apparent to those skilled in the software art. 

10 

The present invention preferably includes a computer program product, 
which includes a storage medium comprising instructions that can be used to direct 
a computer to perform processes in accordance with the invention. The storage 
medium preferably includes, but is not limited to, any type of disk including floppy 
15 disks, optical data carriers, compact discs (CD), digital video discs (DVD), 

magneto-optical disks, read only memory (ROM), random access memory (RAM), 
electically programmable read only memory (EPROM), electrically eraseable 
programmable read only memory (EEPROM), magnetic or optical cards, or any 
type of media suitable for storing information. 

20 

Stored on any one of the above described storage media, the present 
invention preferably includes programming for controlling both the hardware of 
the computer and enabling the computer to interact with a human user. Such 
programming may include, but is not limited to, software for implementation of 
25 device drivers, operating systems, and user applications. Such storage media 
preferabty further includes programming or software instructions to direct the 
general purpose computer to perform tasks in accordance with the present 
invention. 

30 The programming of the computer preferably includes software for 

digitizing and storing images obtained from the image acquisition device (helical 
computed tomography scaimer). Altematively, it should be understood that the 
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present invention may also be implemented to process digital data derived from 
images obtained by other means, such as x-rays and magnetic resonance imaging 
(MRI), positron emission tomography (PET), ultrasound, optical tomography, and 
electrical impedance tomography. 

5 

The invention may also be implemented by the preparation of application 
specific integrated circuits (ASIC), field programmable gate arrays (FPGA), or by 
interconnecting the appropriate component devices, circuits, or modules, as will be 
apparent to those skilled in the art. 

10 

A. LUNG SEGMENTATION 

This section discusses a lung segmentation algorithm in detail. First the 
model of the lungs will be presented with a discussion of the difficulties with this 
1 5 task. Then each stage of the algorithm is presented. 

Referring now to Figure 1, the object of the algorithm is illustrated. The 
top 6 images in Figure 1 are selected slices from a single CT scan, and the bottom 
6 images are the lung regions extracted from those slices. The lungs are modeled 
20 as a low density region surrounded by a high density one. Image filtering is 

preferably first performed to minimize the effects of image noise. A simple linear 
discriminant ftinction based on CT voxel values (photon density) is preferably 
used to partition the scan into lung parenchyma and solid structures. Further 
filtering can be used to compensate for imperfections in the thresholding. 

25 

The photon density of the structures within the lungs may be found by 
manually segmenting the lung into different regions and analyzing the results. 
Referring to Figure 3, a table presents the range and mean values for the intensity 
of solid tissue and lung parenchyma. These values were obtained by manually 
30 selecting regions within a lung CT scan. Regions of at least 100,000 voxels were 
selected from four 2.5mm full lung scans. The distribution of intensities is shown 
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in Figure 2. Note that there is an overlap between the intensity of the lung region 
and the intensity of the surrounding solid tissue. 

Referring again to the table in Figure 3, it is almost possible to partition 
5 lung parenchyma and bone using a linear discriminant function. All voxels with a 
gray level value less than 500 can be considered parenchyma because the 
minimum intensity for soUd tissue is 545. Unfortunately, vessels, nodules, and 
partial voxels can not be partitioned so distinctly. Using the linear discriminant 
function as described will result in classification errors where the histograms in 
10 Figure 2 overlap, most prominantly between gray levels of 550 to 800. All values 
of gray level and intensity in the application are in GE \mits which equates to 
Hounsfield units HU as follows: HU = GE units - 1024, for example 1000 GE = - 
24 HU. 

1 5 One major challenge is handling the streaking artifacts found in 2.5mm 

scans with a low radiation dose and a high pitch, such as those done in screening 
studies. When an X-ray beam passes through a region with a high attenuation 
coefficient, the remaining beam contains a greater proportion of high energy. This 
makes the path-integral of attenuation a non-linear function of distance (J. Hsieh, 

20 "Adaptive Trimmed Mean Filter for CT Imaging," SPIE, Vol 2299, 1994; H. 

Soltanian-Zadeh, J.P Windham and J. Soltanianzadeh, "CT Artifact Correction: 
An Image Processing Approach," SPIE, Vol 2710, 1996 which is incorporated 
herein by reference). In the case of lung CT scans, this artifact is prevalent through 
the shoulders, resulting in noise in the first 25% of the slices. The noise is 

25 characterized as both high intensity and low intensity horizontal streaks running 
through the image. 

Finally, there are other structures which confound the segmentation 
scheme. Air or gas found outside of the lungs will have the same characteristics as 
30 the lung themselves. Specifically they consist of a low density region surrounded 
by a high density one. It is necessary to discriminate between the lungs and other 



21 



wo 03/024184 



PCT/US02/29366 



stray tissue or air pockets. An example of such a region is shown in Figure 4, 
which occurs well below the lungs and the diaphragm. 

The segmentation algorithm consists of creating a mask representing the 
5 lung volume and applying an AND operation between the mask and the original 
image. The algorithm for creating the mask is set forth in Figure 5, and Figure 6 
outlines these steps graphically while the motivation behind this routine is 
described in detail in the following sections. 

10 The main object of image filtering is to remove the streaking artifacts on 

the first 25% of the slices due to x-ray beam hardening through the shoulders. The 
streaking artifacts are characteristically horizontal and tend to greatly distort the 
segmentation. Many of the later steps in this algorithm operate on binary images 
and are quite susceptible to noise. A combination of both median and mean 

1 5 filtering reduces the noise significantly. Median filtering preserves the edges 

around the chest cavity while eliminating much of the streaking effect. Since the 
streaks are horizontal, a small vertical mean filter will blur the streaking regions, 
further reducing their effect. This filtering stage is used only to generate a clean 
image for the mask as the original pixel values are retained in the final segmented 

20 lung image. 

Experiments have shown that a median filter size of 4x4 produce the best 
results. Smaller filters did not reliably eliminate noise caused by the streaking 
artifacts. Further experiments have shown that a mean filter size of 1x3 produces 
25 excellent results. The results of filtering are shown in Figure 9, but the importance 
of this step is seen more clearly in Figure 7 after the thresholding operation. The 
filtering is preferably performed on 2 dimensional slices rather than 3D volumes. 

On large data sets, such as Ml chest CT scans, median filtering is a time 
30 consuming operation. To speed up the algorithm, it is important to note that the 
streaking artifacts appears only at the beginning of the scans. Preferably the 
median filter for the first 25% of the slices is only computed to reduce the 
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computation time of this stage by 75%. The mean filter is much smaller (1x3) and 
the computational gains of only running it on part of the image are insignificant 

Most of the segmentation is performed by a simple threshold operation. As 
5 seen in Figure 2, the distribution of intensities of lung parenchyma and solid tissue 
is bimodal. A simple threshold of 500 can therefore be applied to the entire image. 
For the most part, this will separate the image into 2 regions: (1) solid tissue and 
(2) everything else. Included in "everything else" are the lungs, trachea, air 
surrounding the body, and black background outside of the field of view. By 
10 filling in all pixels which touch the image border and inverting the image, the lungs 
are clearly segmented from the chest cavity, as shown in Figure 8. 

As seen in Figure 2 and the table in Figure 3, there is some overlap between 
the intensities of structures in the lungs and the intensities of solid tissue outside 
15 the lungs. Specifically, there is significant overlap in the range between gray 
levels of 550 and 800. It is therefore impossible to threshold out the ribcage and 
keep all the vessels, as seen in Figure 8. 

Figure 9 illustrates the effect of thresholding both with and without image 
20 filtering. The noisy image on the left shows the result of thresholding without 
filtering and the segmented image on the right shows the result of thresholding 
after filtering. 

Up until this stage in the algorithm, any regions with similar characteristics 
25 to the lungs will be preserved. This includes air and gas found v^thin the body. 
This also includes random noise in the background (i.e., outside the thorax) which 
was not removed by image filtering and thresholding. Selecting only the largest 
component(s) in the image, we will eliminate minor air pockets and all stray noise 
which is prevalent throughout the images. 

30 

Figure 10 shows a region which is removed by retaining only the largest 

component(s) and discarding all other structures. To perform the selection, 
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preferably all regions with voltimes greater than 1% of the total number of pixels 
are chosen in the image, as explained in S. Hu, E. A. Hoffman, and J . M. 
Reinhardt, "Automatic Lung Segmentation for Accurate Quantitation of 
Volumetric X-Ray CT Images," IEEE Transactions on Medical Imaging, Vol 20, 
5 No 6, June 200 1 , which is incorporated herein by reference. 

As previously mentioned, the thresholding operation will segment out 
blood vessels along with the ribcage, resulting in holes and splotches in the lung 
region, as seen in Figure 8, These can be filled in using a morphological filtering 
10 operation. By performing a closing, which consists of a dilation followed by an 
erosion, the spott}^ regions and holes within the lungs can be removed. Some of 
the vessels and other structures are quite large, necessitating the use of a large 
closing kernel. The results of morphological filtering are shown in Figure 11. 

15 Morphological filtering with large structuring elements is computationally 

expensive. Fortunately, large structures that require a large filter kernel only occur 
in the middle third of the slices of the lungs. Therefore, the lungs can be divided 
axially into three distinct regions 1, 2, and 3 as shown in Figure 12. 

20 As the slices corresponding to each region show, preferably a small 

structuring element is used in regions 1 and 3 where only small vessels occur, and 
a larger one is used in region 2, where large vessels are found. Specifically, a large 
2D circular filter is preferably used in region 2 (typically 28 pixels in diameter) 
and filter of half the diameter is used in regions 1 and 3. Referring now to Figure 

25 13, a table that shows timing trials using a varying kemel size versus using a fixed 
large kemel for the entire image. In this experiment, a circular filter of diameter 28 
was used as the fixed kemel size. This was compared to using the same kemel in 
region two and reducing it to a diameter of 14 in regions 1 and 3. 

30 Once the binary mask of the lungs is generated, a simple AND operation 

with the original CT scan will yield the segmented lungs. Finally, the image is 

clipped such that the new bounding box is just large enough to contain the lungs. 
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B, LUNG VOLUME MEASUREMENTS 

This section discusses the algorithms used for lung volume measurements. 
First the motivation behind measuring the lung volume will be explained. Next, a 
5 model of the lung will be presented, followed by the algorithm for calculating the 
volume. 

Once the lung has been properly segmented, the volume of both tissue and 
air in the lungs can be calculated. The objectives of this research are as follows: 

10 

1 . Accurately measure lung volume. 

2. Determine both the inter-patient and intra-patient variation in lung volume. 

3. Explore the effects of inspiration on CAD algorithms. 

4. Compensate for any effects discovered in (3). 

15 5. Measure the relationship between cross-sectional area of the lungs and 

entire lung volume. 

Once statistics of the lung volumes have been compiled the effect of 
inspiration on nodule volume measurements will be studied, A positive correlation 

20 between lung volume and nodule size in a large set of repeat scan will determine if 
such a relationship exists. It is also important to note that full chest scans are not 
always performed during a repeat examination. Typically, a high resolution scan 
of the suspicious region is performed. Such scans contain only a few cross- 
sectional slices of the region of interest. Therefore, the relationship between 2D 

25 cross-sectional area and inspiration must be explored. 

The lungs can be modeled as containing air with a gray level value of zero, 
and tissue with a gray level value of 1000. All intensity values in between that 
range must be a combination of air and tissue, confoimded by the partial voxel 
30 effect. Based on this assumption, a gray level of 500 represents a volume 

comprised of 1/2 air and 1/2 tissue. Equation 1 depicts this formally, where Vis 
volume and /is intensity. 
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•^tissue 



For example, a gray level of 250 represents 25% tissue and 75% air. If the 
image resolution is .5mm x .5mm x 2.5mm per voxel, the volume of a voxel is 

3 3 

.625mm and the volume of air in that voxel is . 1 56mm . 



Referring now to Figure 14, a matrix is generated from the segmented lung 
image in which each column represents a gray level value and each row represents 
10 a slice in the scan to perform this analysis. The entries in the matrix are the 

number of times that a particular gray level appears in that slice of the image. This 
data structure is simply a histogram of the image, separated into slices. With such 
a structure in place, it is easy to perform and repeat the volume measurements on 
both the full lung volume and on an individual slices. 

15 

To calculate the volume in the lungs, the number of voxels in the 
segmented image is counted and multiplied by the volume of each voxel. This 
represents the volume of the entire lung. This can be performed quickly using the 
histogram because the counting has already been completed. Using the model 
20 described above, that each voxel is part air and part tissue, the sum of "tissue 
voxels" multiplied by the volume of a voxel equals the amount of tissue in the 
lung. The equations for calculating the volume of the lungs and the volume of the 
tissue in the lungs are as follows: 



(3) 
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where K represents volume, /represents voxel intensity, and is tiie number 

of pixels of intensity i in slice s of image L. 

Often times with repeat scans, only a high resolution scan of the region of 
5 interest is acquired. As a result, full lung volumes are not always available to 

measure the amount of inspiration in repeat scans. As a solution, the area of single 
axial CT images from two scans may be used to estimate the variation in 
inspiration between the two. The histogram data structure described above can be 
used to calculate the cross-sectional area of a single slice. Given the volume of the 
10 2 entire lungs, it may be possible to relate the cross-sectional area of each to the 
volume of the entire scan. 

The lung segmentation algorithm was run on 196 individual chest scans and 
24 pairs of repeat scans from a Comell University ELCAP database. Out of the 
15 196 scans, 115 contained 2.5mm slices and 81 contained 5mm slices. The 24 pairs 
of repeat scans consisted of 19 5mm scans and 5 2.5mm scans. Repeat scans of 
different reconstructions were not considered. In-plane resolution varied between 
.54mm/pixel and .79mm/pixel. 

20 Images were acquired on a either a General Electric® LightSpeed™ or 

HiSpeed™ Helical CT scanner. X-ray tube current ranged from 40-200mA and 
the tube potential was either 120 or 140 kVp. Testing is preferably performed on a 
dual-processor 700MHz computer. The algorithm ran in approximately 6-8 
minutes on a full lung scan with 2.5mm contiguous slices. Without using the 

25 varying filter sizes as described above, the algorithm ran in about 14-15 minutes. 

C. FINDING LOCATION AND SIZE OF PULMONARY NODULE 

Given an initial seed point inside the nodule, an algorithm finds the center 
30 location of the nodule and the radius of the sphere circumscribing the nodule. The 
following subsections describe the algorithm and choices of functions. 
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The algoriUun is provided a candidate location that is some seed point P, 
which is within the nodule in the image, and a starting radius r and delivers the 
optimal location F and radius of the nodule. The algorithm for finding the 
location P' and radius r' of the nodule is shown in Figure 15. The value of r is 
5 incrementally increased and for each value of r the best location of the nodule 
center P is computed, A new location is calculated by maximizing the correlation 
(or the response) between the image and a locater template function. The template 
function is chosen so that the maximum response is given when the template is 
centered over a nodule with radius r as discussed further below. Thus, maximizing 
10 the correlation yields the most likely location of a nodule with the current radius. 

The algorithm determines the radius of the nodule by terminating itself at 
the appropriate iteration. Another template function, centered at P, is applied to 
the image at the end of each iteration. The sizing template function is designed so 
1 5 that its value decreases dramatically when the radius of the function has become 
larger than the radius of the actual nodule. Thus, when there is a large decrease in 
the value of the second template function, the radius of the actual nodule has been 
exceeded and the algorithm terminates and returns the location and radius from the 
previous iteration. 

20 

The nodule finding algorithm was developed assuming a two level model 
where the nodule, blood vessels, and pleural wall are all a single high intensity, 
while everything else is a single low intensity. The template functions use 
parameters of location and radius. The radius does not refer directly to the size of 
25 the template function. Rather, the radius is a parameter that describes the size of 
the nodule for which the template will detect; meaning that the template does not 
necessarily need to be zero past the extent of the radius 



The locater template functions is designed to locate the approximate center 

30 of mass of a nodule (for nodules attached to the pleural surface a repeatable central 

location is obtained since the location of the actual center of mass is not 

knowable). The function has four parameters: the center in x, y, and z, and the 
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radius. In the algoritibm, the radius is fixed at each iteration, but the center of the 
template can move. The locater template function was designed, with a given 
radius, so that it is maximized when the function is centered over a nodule of the 
same radius. 

5 

Consider a 1-D example where an isolated nodule is modeled as a large 
region of high intensity. Figure 16A shows an example model where the nodule is 
the center shaded region, while the outer regions consist of noise (composed of 
high-intensity objects like vessels). 

10 

Consider some simple 1-D template functions shown in Figure 17. The 
square template function. Figure 17 A, biases everything equally, causing the 
template to be very sensitive to noise in the periphery. The sensitivity to noise can 
be reduced by tapering the template function to be shaped like a triangle, Figure 

15 17B, or a Gaussian, Figure 17C. The large values near the center will keep the 
function centered over the high-intensity region of the nodule, while the lower 
values further from the center will diminish the effect of any other high-intensity 
objects. Figures 18A through 18C illustrate where the Gaussian template function 
will be maximized for cases where the size of the function is smaller, the same 

20 size, and larger than the size of the nodule. 

Now consider a 1-D model of a nodule attached to the pleural wall, shown 
in Figure 18B. The nodule is the lightly shaded central region, the pleural wall is 
the darkly shaded region on the right, and the lightly shaded region on the left is 
25 some noise (blood vessel). Although shaded differently, note that the nodule and 
the pleural wall are at the same intensity level. With this model, the Gaussian 
template function will fail to perform correctly. The template does not prevent 
itself from slipping into the wall and centering itself over the combination of the 
nodule and the wall. 

30 

Nodules on the pleural wall can be located by choosing a function that 
forces itself away from Hie wall. This can be achieved by using a function that has 
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positive values within the radius of the functioii and negative values outside the 
radius, like Figures 17D through 17F. Two examples of this type of function are 
the Laplacian of the Gaussian (LOG) and the difference of Gaussians (DOG). 
When maximizing the response of the template function, there are two forces 
5 fighting against each other. The positive portion of the template function will try 
to keep the function centered over the high-intensity regions. Countering this, the 
negative portion of the template will try to push the template away from the high- 
intensity regions, resulting in. an equilibrium along the edge of the high-intensity 
region of the nodule. Figures 18D through 18F show how the LOG-like template 
10 function should behave on the 1-D model of a nodule on the wall. If the total 
weights, the distribution, and the sizes of the positive and negative regions are 
balanced correctly, then the LOG-like template function should behave 
accordingly. 

1 5 The LOG-like template function can also be used on isolated nodules. The 

positive region of the template function will stay near the center of the nodule, 
while the negative regions will force the function away from any noise in the 
periphery of the image. Figures 18G through 181 show how the LOG-like function 
will behave when applied to the ID model of an isolated nodule. When the size of 

20 the template is smaller than the nodule, the template will hug the edge of the 
nodule and when the positive region of the template is larger than the nodule, it 
will be centered over the nodule. 

The above illustrates the preferred behavior of the locator template function 
25 for the ID model of an isolated nodule and a nodule on the wall. The model is 
now extended into two and three dimensions. For the 2D case, the ID LOG-like 
function is extended into two dimensions. The template now consists of two 
circular regions; the smaller region contains positive values and the region outside 
the smaller region contains negative values. Figure 19 shows a few examples of 
30 how the 2D template function is expected to behave for isolated and pleural 

nodules. For the 3D case, the locater template function will have spherical regions 
with an analogous distribution as the 2D and ID cases. 
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One thing to note about the locater template function is that it does not need 
to have a constant total weight for different radii. The algorithm does not 
maximize the locater template function throughout all the iterations and all the 
5 radii; instead, it is maximizing the template function for a fixed radius in each 

iteration. Thus, the template is not required to have a constant weight for different 
radii. It should suffice to only maintain the proper scaling between the sizes and 
the distribution of the template regions. 

10 The sizing template function is used to indicate when the radius of the 

function is equal to the size of the sphere circumscribing the real nodule. The 
sizing template function uses the radius during the current iteration and the 
location of the candidate nodule, determined by the locater template function. 
After each iteration in the algorithm, a value is calculated using the sizing template 

15 function. Two preferred methods for determining when the algorithm has found 
the location and size of the nodule in the image are discussed below. 

The first method is to look for a large change in the response of the sizing 
template function. Using the LOG-like locater template function, the positive 

20 region, the inner circle, is expected to be overlayed on the high-intensity values of 
the nodule. As the radius of the template function increases, the positive region 
will remain within the nodule. This is illustrated in Figure 20. Eventually the 
positive region will become larger than the nodule and there will be a sharp 
decrease in the response of the sizing template function between iterations of the 

25 algorithm. To find this stopping criterion, a sizing template function is chosen so 
to be a mean filter shaped like a sphere. When the sphere becomes larger than the 
nodule, the distribution of values inside the sphere changes, causing the response 
of the template function to change dramatically. 

30 Figure 22 shows the response of the mean filter fimction to various sizes of 

ideal nodules (spheres of intensity 1000). For the ideal nodule, the response of the 
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filter function is constant until the radius of the filter fiinction exceeds the radius of 
the ideal nodule. When this happens, the response of the filter quickly decreases. 

The second method of finding the size of a nodule is to design the sizing 
5 template function so that it has a specific response when the template function is 
the same size as the real nodule. For example, if a 3D Gaussian function, where 
the standard deviation is equal to the radius, has the response shown in Figure 23. 
When the response of the Gaussian template is 200, the radius of the Gaussian is 
the same as the radius of the nodule, thus providing a stopping criterion for our 
10 algorithm. 

One very important property of the sizing template fiinction is that it must 
have a constant total weight regardless of the size of the template. This property 
will insure that the response of the template will be consistent between different 
1 5 sizes of nodules. 

The implementation offhe algorithm will now be discussed. This will 
include a more detailed description of the algorithm, a mathematical definition of 
some template fimctions, calculation of the initial seed point, introduction of two 
20 search methods for the maximization of the locator template function, quick 
function evaluations, and the termination criterion. 



A detailed description of the algorithm is shown ia Figure 21 . The 
algorithm is modified so that it terminates when the response of the sizing template 
25 function is within e of the target value. This feature is added to the algorithm so 
that the radius does not go too far past the actual radius of the nodule. 

The algorithm takes in a 3d CT image of a nodule. First, the image is 
clipped at intensities over 1000. This is done to minimize any affect that the very- 
30 high intensitied bone will have on the algorithm. 
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The next step is to define the locator and sizing template functions. The 
template functions are preferably either the Gaussian, the Laplacian of the 
Gaussian, or the difference of Gaussians. The template function has four 
parameters, (sx, sy^ sz^ a), which correspond to the x-location, y-location, z- 
5 location, and radius. The initial location is either calculated from the image or 
specified by the user, and dependiag on the type of sizing template function, liie 
target response value is set. 

Starting with the initial conditions, a radius of 1, and a radius step-size of 1, 
10 the algorithm finds the best locator template response by using a search method to 
move the center of the template function. Preferably the search method is a hill- 
climbing method as discussed below. Once the maximum value is found for the 
given radius, the radius is iacreased by the radius step-size and the search is 
repeated using the current location as the starting point. 

15 

At each iteration the response of the sizing template function, y is 
calculated. The algorithm iterates until the sizing response is less than the target 
value. If Y is within 6 of the target, then the algorithm is finished. Otherwise, the 
algorithm backtracks by setting the location and radius parameters to their values 
20 in the previous iteration, dividing the radius step by two, and repeating the search 
procedure. Finally, after the target value is reached within e, or the radius step-size 
is smaller than some a, the search process terminates and the current parameters 
are the location and radius of the nodule. 

25 The template function is a function of four parameters: the center 

(location) of the template, {sx^sy^sz)^ and the radius of the template, a. The response 
of the template function is calculated by taking the correlation between the image 
im and the template M as in Equation (4). 



y vy 

30 T{sx, sy, sz, a) = zm(x, y, z) * M^,,^,^ (x, y, z) (4) 
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The locater template function is designed to produce the largest response 
when the kemal is centered over the target object. Thus, by maximizing this 
function at each iteration of the nodule finding algorithm, the best candidate for a 
nodule of radius r is found. There are many types of functions to choose iBrom, but 
the preferred functions are the Gaussian, the Laplacian of the Gaussian, and the 
Difference of Gaussiaus. 



The 3D gaussian is a strictly positive symmetric function defined by 
Equation (5). The Gaussian function can be used to find nodules when the pleural 
surface is not present in the image. The metric function will have its highest value 
when it is centered over the nodule. A graph of the Gaussian is shown in Figure 
24. 



The Laplacian of the Gaussian (LOG) is defined in Equation 7. This 
function has a positive weight close to the center and a negative weight further out, 
as seen in Figure 24. This function is useful in locating nodules on or near the 
pleural surface. The positive interior will latch onto the nodule, while the negative 
exterior will keep the kemal from moving towards the wall. In the algorithm, the 
LOG has the disadvantage that the initial location must be within the nodule. 
Otherwise, the negative weight of the exterior will cause the LOG to grow away 
from the nodule. 



LOG, 



^sx,sy,sz,a{x,y,z) 



■ + - 



dx^ dy^ dz' 



(6) 



3-4r-^-4T I * (7) 
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The Difference of Gaussians, given in Equation (8), is similar to the 
Laplacian of the Gaussian in that they both have positive and negative regions. 
However, changing the o"^, a^, m^, and values will lead to different sizes and 
weights of the two regions, allowing for a more customizable template than the 
5 LOG. 



The initial seed point should be somewhere within the nodule in the image. 
10 This may be provided by user input, or it can be estimated from the image. A 

simple way to calculate the initial location would be to use the center of the image: 



center^ = im.xlo-\r(imjchi-im.xlo+l)/2 
center-y = im.ylo-^(im.yhi'im.ylo+iy2 
15 center^ = im^lo+Qm^hi-irn.zlo'^iyi 



When using the Gaussian for the locator template function this may be fine, 
but when using the LOG or DOG functions the initial starting location must be 
within the nodule. For these functions, the center of mass (COM) of the image is 
20 preferred for setting the initial conditions. The COM is calculated by first 

thresholding the image at half of the largest pixel value in the image and then 
calculating the center of mass using the standard equations: 



25 

COM _ ^ ^ ^ 

' Z S Z uj.M'^J^^"> 
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ZEE 



im{ij,k) 



The method used to search for the optimal nodule location has a major 
impact on the running time of the algorithm. Any kind of fancy search procedure 
could be used, but this is not necessary for this application. From iteration to 
iteration the new optimal location will not deviate much from the previous optimal 
location. Since the previous optimal location is used as the starting point for each 
search, the search method will not have to look very far for the maximum in the 
search-space. 

Fancy search methods, such as Powell's method described in W. Press, 
Numerical Recipes in C, 2nd Edition, Cambridge University Press, 1992, which is 
incorporated herein by reference, are usually less efficient when moving small 
distances because these algoritlmas are designed to search a very large space. In 
other words, they are optimized to tradeoff between the efficiency in moving giant 
steps and the efficiency of moving in small steps. 

In the end, hill-climbing, a greedy search method, is preferred for the 
nodule finding algorithm because it is more efficient when the optimal location 
doesn't move very far. The Hill-Climbing is shown in Figure 25. Using an the 
initial location and a stepsize, the hill-climbing algorithm evaluates the locator 
template function for each of the possible 6 moves in parameter-space and then 
moves in the direction of the largest decrease. The algorithm terminates when there 
is no move that decreases the metric. 

The search procedure requires many evaluations of the correlation between 
locator template function and the image. These calculations can be costly because 
the template function usually involve exponentials. Fortunately, in any given 
iteration of the search procedure, the radius of the template function is held 
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constant. It is possible to save calculations by storing the template function values 
in a temporary image (the classic engineering dichotomy of time versus space). 

The values of the locater template function centered at the origin are stored 
5 in a temporary image at the beginning of each search procedure. For the Gaussian 
and the LOG, values are only calculated out to 3 or 5 times the radius a because 
the functions converge to zero. Also, the function is only evaluated at positive 
location points because these functions are radially symmetric, 

10 Now when the response of the template function is required, the temporary 

image is offset by the location of the template function, and the correlation 
between the image and the temporary image is calculated. 

The target response determines when the algorithm should terminate. The 
1 5 target value is chosen such that the algorithm will terminate when the location and 
radius parameters best circumscribes the nodule. The target value can be found 
either through an ideal model of a nodule or through experimentation. Figure 23 
shows the response of a Gaussian function to several ideal nodules (spheres) of 
different radii. In the ideal case, the target response of the Gaussian sizing 
20 function should be 200. The best target response has been determined to be 300 by 
experimentation with a small subset of real nodules from the database. Using a 
small number (four) of nodules, the response of the sizing function was tracked 
and the value that causes all the nodules to be circumscribed was selected. 

25 Based upon the experimentation, the locator template functions shown in 

Figure 26 will be stable. The large negative values near the border between the 
positive and negative regions help the template find the edges of the nodule. In 
addition, the difference of Gaussians (DOG) can be used to correctly balance the 
weights and sizes of the positive and negative regions in the template so that the 

30 template will not slip into the pleural wall 
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The nodule finding algorithm is preferably implemented in the C 
programming language for a VisionX software package on a FreeBSD system. 
VisionX software provides computer tools and programs for the analysis and 
visualization of image data. It is suitable for a wide range of image analysis 
5 applications and is designed to address the processing needs of multidimensional 
image sets that arise both from temporal image sequences and from image 
modalities that involve three-dimensional data collection. 

VisionX software has been used in a wide range of research applications 
10 including multispectral image analysis, three-dimensional object recognition, 
multiframe image analysis, target tracking, neural networks, biological cell 
analysis, and three-dimensional biomedical image analysis. Important features of 
the VisionX software include the ability to handle multidimensional image sets, a 
wide range of available processing fimctions, and a flexible tagged data format that 
1 5 facilities the automatic recording of tihie history of a file. 

FreeBSD is an advanced operating system for x86 compatible, DEC Alpha, 
and PC-98 architectures. It is derived from BSD UNIX, which is a version of 
UNIX developed at the University of California, Berkeley. FreeBSD offers 
advanced networking, performance, security, .and compatibility features. FreeBSD 
20 can be installed from a variety of media including CD-ROM, DVD-ROM, floppy 
disk, magnetic tape, an MS-DOS partition, or if there is a network connection, it 
can be installed directly over anonymous FTP or NFS. 

The nodule finding algorithm preferably inputs a floating point, byte, or 
25 short image and outputs the location and size of the nodule. Many options are 
available, including translation and radius search limits, a choice between search 
strategies (Hill-Climbing or Powell's method), initial starting location, type of 
template ftinction, and termination criteria. The program also preferably outputs 
several types of images for debugging purposes. The program was expanded to 
30 allow the algorithm to run on two-dimensional images. 
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D, PULMONARY NODULE REGISTRATION 



The algorithm registers puhnonary nodules in 3-d images of high-resolution 
focused CT-scans. For the most accurate characterization it is important to have 
5 images from high-resolution CT scans because these images will better localize the 
boundaries of the nodule. A rigid-body transformation model is assumed, meaning 
that in general the structures in the images are confined to simple translation and 
rotation. In order to simpUfy the transformation model, an isotropic image-space is 
also assumed. Using these assumptions, an algorithm was developed to register 
10 pulmonary nodules from two different time periods by defining a metric between 
the two images and performing a minimizing search on the metric. 



The algorithm requires two input images, and produces an output that is the 
first image registered to the second image. The rigid-body registration algorithm is 

1 5 shown in Figure 27. The two input images, im^ and im^, are converted to floating- 
point pixel format. Preferably the images are masked to ignore irrelevant pixel 
data (e.g. bones). Next, the initial conditions for the rigid-body transformation is 
determined. The metric between two images is defined as a number that represents 
how closely two images are related to each other. By selecting the appropriate 

20 rigid-body parameters, the metric between the second image and the 

transformation of first image can be minimized or maximized depending upon the 
registration metric, resulting in the best registration for the given metric. 



The first step in the registration algorithm is to perform some pre- 
25 processing on the input images. If necessary, the input image is converted from 
byte or short pixel-formats to the floating point pixel-format. 

If mask images are provided, the input images are masked by setting the 
appropriate pixels to the background value (-1000). Masking an image tells the 
30 registration algorithm to ignore certain areas of either image when attempting to 
perform the registration. This can be usefiil if a large structure in an image is 
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causing misregistration of the object of interest. By using a mask, the algorithm 
can be told to ignore the larger structure, thus allowing the object of interest to be 
registered to itself. 



5 The registration metric is defined as a function of two images, un^ and im^ 

and six rigid-body transformation parameters, {t^^ty^t^^r-^yTy^r^. To calculate the 
metric, the first image is transformed using the rigid-body parameters, resulting in 
im^^. One definition of the metric value is the correlation between the transformed 
image im^^ and the second image ini^^ as given in the following equation. 

10 

^^=T7 ZEE imxtihm^ iimiSJM) (9) 

imu{ij\k) 9^ackground 
i^iQj^k) j^ackground 

15 where i\r is the number of pixels that are not background pixels in either un^^ or im^. 



The correlation metric is not absolute in the sense that there is no absolute 
best correlation that indicates that the two images are registered perfectly. The 
correlation is dependant on the distribution of the pixels; images with larger pixel 
20 values will produce a larger correlation. An absolute metric allows for a better 
sense of how well two images are registered to each other regardless of pixel 
distributions. 

Instead of correlation, a mean-squared-difference (MSEJ) metric, defined in 
25 the following equation, can be used: 

= ^ E ZZ (i^uihj.k) - im2iij\k)f (10) 

imu(ij\k) T^ackground 
im2(ij,k) background 
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10 



where iV is the number of pixels that are not background pixels in either im^^ or im^. 
With the MSD metric, perfectly registered images will produce a metric of zero. 

The rigid-body transformation uses six parameters: (Jx,ty,tz,rx,ry,rz); 
translation in x, translation in y, translation in z, rotation about the x-axis, rotation 
about the y-axis, and rotation about the z-axis. The transformation is a mapping of 
a point V in 3-d space to a point v' in transformed space defined by the following 
equation: 



where R^^ Ry, and R^ are the rotation matrices defined as: 



(11) 



15 



1 0 

0 cos(rJ 
0 sin(r,) 

cos(ry) 
0 

-sin(r^) 

cos(r^) 
sin(rj 
0 



0 

-sin(rj 
cos(r,) 

0 sia(ry) 

1 0 
0 cos(r^) 

-sm(rj 0 
cos(r,) 0 
0 1 



(12) 



(13) 



(14) 



The first image is transformed so that it has the same boxmding box as the 
second image. The transformation uses eillier linear interpolation or nearest- 
20 neighbor interpolation, and pixels that are transformed firom outside the bounding 
box of the first image are set to the background value. 
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The initial rigid-body transformation parameters must be determined before 
a minimization search can be performed. The initial rotation parameters are all set 
to zero because only a small amount of rotation is expected between the two 
images. The initial translation parameters should be set so that the two nodules 
5 will overlap at the first iteration of the search procedure. Assuming that the 

nodules are in the center of the images, the initial translation parameters can be set 
to the difference between the centers of the two images, where the center of an 
image is defined as: 



10 center-^ = imjclo'^(im,xhi-im.xlo+l)/2 

centery = im .ylo+(im.yhi-im.ylo'^iy2 
center^ = Un,zlo^{pn,zhi''imjzlo^\yTL 



This assumption is not always true and a better method of determining the 
15 initial translation parameters is to use the centers of mass of the two nodules. Both 
images are thresholded at half of the largest pixel value in the image and the center 
of mass is calculated using the standard equations: 



r^n^^ ^ II ,,jjm(i,j,k)^i 
^O^x = ^ ^ ...... 

20 COM ^ ^ 



ZEE t.j.M''J'''^ 
Z Z Z ,jM''J'^^*^ 

COM = '-^^ 

' Z Z Z .,jJrn(iJ,J^ . 



The initial translation parameter is equal to the difference between the centers of 
mass of the two images. 

25 

Using the registration metric, a function of six parameters, and the initial 
conditions, the registration problem is now converted into a 
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minimization/maximization search problem. Specifically, given some initial 
condition, a minimum or maximum of the metric function is found by searching 
- the 6-dnnensional parameter space. The following subsections describe three 
different search methods: exhaustive search, hill-climbing, and Powell's method. 

5 

The exhaustive search is a brute force method that finds the minimum (or 
maximum) by calculating the metric for every possible set of parameters. Of 
course the number of possible parameters in a 6-dimensional search space is 
infinite. In order to reduce the search space, a limit of translation and rotation 
10 deviation from the initial condition, along with rotation and translation stepsizes, is 
imposed. Even with this search space reduction, the exhaustive search is still 
impractical. 

Consider a case where translation is limited to 5 pixels in each direction 
15 and the rotation by 5 degrees in each direction. With stepsizes of 1 pixel of 

6 

translation and 1 degree of rotation, the exhaustive search will perform 1 1 or 
1,771,561 metric evaluations. Suppose each metric evaluation takes 0.1 seconds, 
then the entire algorithm will take around 170,000 seconds, or 2,833 minutes, or 
just under two days. Even with such a seemingly small search space and 
20 appropriate step-sizes, the exhaustive search will take days to finish! 

Hill-Climbing is a greedy search method. Starting at the initial condition 
and a set of translation and rotation stepsizes, the algorithm evaluates the metric 
for each of the possible 12 moves in parameter-space and then moves in the 
25 direction of the largest decrease (or increase). The algorithm terminates when there 
is no move that improves the metric. 

By using a small step-size it is possible to calculate the best parameters up 
to some precision. Unfortunately, as the step-size is decreased, the running-time of 
30 the algorithm increases because the algorithm must take smaller steps, and thus 
perform more metric evaluations to get to the minimum. One solution to this 
problem is to repeat the algorithm several times while using a decreasing stepsize. 
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A detailed description of the complete hill-climbing algorithm is shown in Figure 
28. 



Powell's method is a multi-dimensional direction-set search algorithm. 
Starting with an initial set of 6 directions and an initial condition, in each iteration 
the algorithm minimizes the metric by moving in each of the six directions. For a 
given direction, any line minimization technique could be used, however Brent's 
Method is chosen because it is a parabolic minimization technique that does not 
explicitly use derivatives. At the end of each iteration (one pass through each of 
the six directions), the oldest direction is replaced by the total direction moved 
during the current iteration. By doing this, the algorithm adapts itself to move 
along the most minimizing path. A more in depth description of Powell's method 
can be found in Chapter 10.5 in Numerical Recipes in C, W. Press, 2"^ Ed., 
Cambridge University Press, 1992, which is incorporated herein by reference. 

Powell's method terminates when either all the parameters in one iteration 
change within some epsilon, or the metric after an iteration does not change by 
more than a tolerance value. 

The registration algorithm works well on images with isolated nodules, but 
has trouble with images where the pleural surface is present because the pleural 
surface becomes the largest structure in the image. The largest structure has the 
largest effect on the metric value, causing the minimization technique to try to 
register the pleural surface to itself. Because the pleural wall can move and change 
shape with inspiration or position, there is no guarantee that the nodules near or 
along the pleural wall be registered correctly if the pleural wall is registered. A 
first method for addressing this problem is to use image masking or pixel masking. 
A mask is created that tells the algorithm to ignore the pleural surface and its 
contents will force the algorithm to register only the nodules together. A second 
method involves applying a function on the pixels so that all the pleural features 
become one intensity. This causes the registration algorithm to ignore structures 
like the ribs. A third method involves reducing the search space. A nodule- 
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localization algorithm is used to find the centers of the nodules in both images. 
Taking the difference of the centers results in the translation parameters that 
overlap the nodules. Finally, the orientation of the nodules is determined by using 
the registration algorithm on only the rotation parameters. 

5 

The rigid-body registration algorithm is preferably implemented in the C 
programming language using the VisionX software library for the FreeBSD 
enviroimient The program, vSregrb, preferably inputs two input images of byte, 
short, or floating point pixel-types and a variety of options. The program 
10 preferably outputs the first input image registered to the second input image and a 
difference image associated with the two images. Both images preferably have the 
same bounding box as the second input image. The rigid-body transformation 
parameters for the registration are preferably stored with the metric and timing 
information with a history of the output images. 

15 

In order to increase the throughput of the registration process, the three- 
dimensional rigid-body transformation is preferably written into the program. The 
code for the transformation is preferably used from the existing stand-alone rigid- 
body transformation program (v3regrb). The rigid-body function preferably 
20 supports both linear and nearest-neighbor interpolation. The program may be 
modified to perform three-dimensional translation-only registration, two- 
dimensional rigid-body registration, and two-dimensional translation-only 
registration. 



25 E. REMOVAL OF THE PLEURAL SURFACE FROM 

JUXTAPLEURAL NODULES IN THRESHOLDED HIGH-RESOLUTION 
CT IMAGES 

The algorithm takes a binary (thresholded) image as its input. The 
30 algorithm also requires the location of a point near the center of the nodule. If the 
image is some region-of-interest that has been generated by a radiologist, then it is 
assumed that the nodule is in the center of the image. In this case the initial 



45 



wo 03/024184 



PCT/US02/29366 



starting point will be the center of the image. The pleural-surface removal 
algorithm is shown in Figure 29. 

The pleural surface removal algorithm works by iteratively moving a plane 
5 towards the pleural-surface. The algorithm starts with an initial point P' inside the 
nodule and a direction towards the wall, d\ The direction is calculated by taking 
the difference between the center of mass of a spherical region centered on P* and 
the starting location P'. 

10 With each iteration, a new location Pf is calculated by stepping in direction 

d from the previous location P/.j . Plane normal to direction d and passing 

through point P, separates the nodule from the pleural surface. This cut nodule is 
the connected component region that contains point P and that is behind the plane 
A. As the plane moves towards the pleural wall, the size of the cut nodule 
1 5 increases. Figures 30A and SOB show the cut nodule for two iterations. 

The algorithm keeps track of the difference A in the size of the cut nodule 
between iterations. When the plane intersects the pleural surface, as seen in Figure 
30C, the difference A increases dramatically. This condition is manifested by the 
20 change in the difference increasing by more than y^iax- When this occurs, the size 

of the cut nodule is minimized by reorientating the plane by changing the direction 
d while keeping the point P fixed, as in Figure SOD. Additionally, the size 
difference and change in difference are recalculated. If the change in difference is 
less than yjnax^ then the algorithm keeps iterating. Otherwise, the algorithm 

25 terminates and returns the cut nodule formed by using the plane in the previous 
iteration. 

Figure 30 shows a two-dimensional example of the algorithm. The cut 

nodule is shown as dark gray while the rest of the nodule and the pleural surface 

30 are light gray. The initial starting point and the initial direction are shown in 

Figure SOA. Figure SOB shows an iteration where the plane does not intersect the 
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pleural wall. Figure 30C shows an iteration where the plane intersects the pleural 
surface, causing the change in the nodule size to increase. In Figure SOD the plane 
is reorientated to minimize the cut nodule, forming a new direction. After another 
iteration. Figure 30E shows the plane intersects the pleural wall again, but 
5 reorientating the plane still leads to a large increase in nodule size, as seen in 
Figure 30F. Thus, the algorithm terminates and retums the cut nodule from the 
previous iteration, Figure SOD. 

In this section, a specific implementation of the algorithm is discussed. A 
10 seeded region growing algorithm is used to find the cut-nodule region, and a hill- 
climbing search is used to perform the minimization when reorientating the plane. 
This section first discusses the equations used for the plane, followed by a 
description of the region growing and hill-climbing algorithms. 

15 The representation of the plane has four parameters and is defined as: 



ax + by-^ cz-\-d=0 (15) 

Given a direction normal to the plane, v, and a point on the plane, P, the 
20 parameters of the plane are calculated as: 

a = (16) 

b =^ Vy (17) 

25 c = (18) 

d = -{aP^+bPy+cPz) (19) 
Finally, a point p is behind the plane if: 

30 apjc '^bpy + cpz + d<Q (20) 
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A recursive region growing is used to determine the cut-nodule region. The 
starting point of the region-growing algorithm is the initial starting location P\ and 
the plane A is calculated from Pf and dp A description of the region growing 

algorithm is shown in Figure 3 1 . For the given pixel the algorithm marks it as 
5 visited and then determines if it is behind the plane A, If it is not, then the 

algorithm exits. Otherwise, the pixel p is marked as part of the cut-nodule region, 
and the algorithm is recursively called for each of the six possible 1 -pixel moves 
from p if the new point has not been visited yet and it is also a foreground pixel. 

10 A simple greedy search algorithm is used to reorientate the plane so that it 

minimizes the cut-nodule size. A new plane is calculated by changing the 
direction normal to the plane. The hill-climbing algorithm is shown in Figure 32. 
The hill-climbing algorithm takes the initial direction and a small stepsize. Next, 
the size of the cut-nodule is calculated for the planes formed by moving the 

15 direction in the six possible coordinate directions. If a smaller size is calculated, 
then the direction is moved in the direction of the largest decrease. The new 
direction is normalized to one, and the algorithm is repeated. If there is no 
decrease in the cut-nodule size, then the algorithm terminates. 

20 The algorithm is preferably implemented in the C programming language 

for the VisionX software package on a FreeBSD UNIX system. The program 
preferably operates on binary images that have floating point, byte, or short pixel 
types. The program preferably outputs the segmented nodule and two debugging 
images. The initial conditions, initial direction, and termination criteria are 

25 preferably specified by options. 

Thus, while there have been described what are presently believed to be the 
preferred embodiments of the invention, those skilled in the art will realize that 
changes and modifications may be made thereto without departing from the spirit 
30 of the invention, and is intended to claim all such changes and modifications as fall 
within the true scope of the invention. 
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WHAT IS CLAIMED IS; 

1 . A method for generating a lung mask for segmenting a lung image jfrom 
voxel data containing the lung image, noise, solid components and surrounding 
background information, the method comprising: 

(a) applying a median filter and a mean filter to the voxel data to reduce 

noise; 

(b) thresholding the noise reduced voxel data to identify solid 
components and other structures; 

(c) identifying the surrounding background in the noise reduced voxel 

data; 

(d) deleting the surrounding background from the noise reduced voxel 

data; 

(e) labeling coimected components of the voxel data; 

(f) determining the largest coimected components to select a Ixmg 
region from the voxel data, the lung region having a geometric form; and 

(g) morphological filtering the voxel data to refine the geometric form 
of the lung region. 

2. The method as defined by Claim 1, wherein the voxel data is associated 
with a plurality of slices through a patient's lung, the slices beginning at about the 
patient's shoulders; and 

step (a) is performed for the first 25 percent of the plurality of slices only 
whereby the computation time is substantially reduced. 

3. The method as defined by Claim 1, wherein the median filter has a size of 
4x4. 

4. The method as defined by Claim 1 , wherein the mean filter has a size of 
1x3. 
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5. The method as defined by Claim 1, wherein the voxel data in step (b) is 
thresholded at a gray level of about 500. 

6. The method as defined by Claim 1, wherein the largest connected 
components are associated with more than about 1 percent of the voxel data. 

7. The method as defined by Claim 1, wherein the voxel data is associated 
with a plurality of slices, the slices being divided into a first end region, a middle 
region, and a second end region; and 

the morphological filtering being performed with 2D circular filter having a 
first diameter in the first end region and the second end region; and the 
morphological filtering being performed with 2D circular filter having a second 
diameter which is about twice the first diameter in the middle region. 

8. A method for measuring lung volume firom a segmented lung image 
obtained firom a scan which includes a plurality of slices of voxel data, the voxel 
data having a gray level value and a volume associated therewith, the method 
comprising: 

generating a matrix of entries fi-om the segmented lung image, the matrix 
having a plurality of columns and a plurality of rows, each of the plurality of 
columns representing the gray level value and each of the plurality of rows 
representing one of the plurality of slices in the scan, the entries corresponding to 
the number of times that the gray level occurs in the corresponding slice; 

determining a number of voxels in the segmented lung image from the 
matrix entries; and 

multiplying the number of voxels in the segmented lung image by the 
volume of each voxel. 

9. A method for measuring volume of tissue in a segmented lung image 
obtained firom a scan including a plurality of slices of voxel data, the voxel data 
having a gray level value and a volume, the method comprising: 
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generating a matrix of entries from the segmented lung image, the matrix 
having a plurality of columns and a plurality of rows, each of the plurality of 
columns representing the gray level value and each of the plurality of rows 
representing one of the plurality of slices in the scan, the entries corresponding to 
the number of times that the gray level occurs in the corresponding slice; 

determining a sum of tissue of voxels in the segmented lung image from the 
matrix entries; and 

multiplying the sum of tissue of voxels in the segmented lung image by the 
volume of each voxel. 

10. The method as defined by Claim 9, wherein the sum of tissue of voxels is 
calculated by summing the product of each matrix entry multiplied by the 
corresponding gray level value divided by a gray level value assigned for tissue. 

11. A lung mask generating apparatus for segmenting a lung image from voxel 
data containing the lung image, noise, solid components and surrounding 
background information, the lung mask generating apparatus comprising: 

a masking unit configured to: 

(a) apply a median filter and a mean filter to the voxel data to reduce 

noise; 

(b) threshold the noise reduced voxel data to identify solid components 
and other structures; 

(c) identify the surrounding background in the noise reduced voxel 

(d) delete the surrounding background from the noise reduced voxel 

(e) label connected components of the voxel data; 

(f) determine the largest coimected components to select a lung region 
from the voxel data, the lung region having a geometric form; and 

(g) morphologically filter the voxel data to refine the geometric form of 
the lung region. 
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12. A lung mask generating apparatus as defined by Claim 1 1, wherein the 
voxel data is associated with a plurality of slices through a patient's lung, the slices 
beginning at about the patient's shoulders; and 

step (a) is performed for the first 25 percent of the plurality of sUces only 
whereby the computation time is substantially reduced. 

13 . A lung mask generating apparatus as defined by Claim 1 1 , wherein the 
median filter has a size of 4x4. 

14. A lung mask generating apparatus as defined by Claim 1 1 , wherein the 
mean filter has a size of 1x3. 

15. A lung mask generating apparatus as defined by Claim 1 1 , wherein the 
voxel data in step (b) is thresholded at a gray level of about 500. 

16. A lung mask generating apparatus as defined by Claim 1 1 , wherein the 
largest connected components are associated with more than about 1 percent of the 
voxel data. 

17. A lung mask generating apparatus as defined by Claim 1 1, wherein the 
voxel data is associated with a plurality of slices, the slices being divided into a 
first end region, a middle region, and a second end region; and 

the morphological filtering being performed with 2D circular filter having a 
first diameter in the first end region and the second end region; and the 
morphological filtering being performed with 2D circular filter having a second 
diameter which is about twice the first diameter in the middle region. 

18. An apparatus for measuring lung volume fi*om a segmented lung image 
obtained fi:om a scan which includes a plurality of slices of voxel data, the voxel 
data having a gray level value and a volume associated therewith, the apparatus 
comprising: 
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a Ixing volume measuring unit configured to: 

generate a matrix of entries from the segmented lung image, the matrix 
having a plurality of columns and a plurality of rows, each of the plurality of 
colxmms representing the gray level value and each of the plurality of rows 
representing one of the plurality of slices in the scan, the entries corresponding to 
the number of times that the gray level occurs in the corresponding slice; 

determine a number of voxels in the segmented lung image from the matrix 
entries; and 

multiply the number of voxels in the segmented lung image by the volume 
of each voxel. 

19. An apparatus for measuring volume of tissue in a segmented lung image 
obtained from a scan including a plurality of slices of voxel data, the voxel data 
having a gray level value and a volume, the apparatus comprising: 

a lung tissue measuring unit configured to: 

generate a matrix of entries from the segmented lung image, the matrix 
having a plurality of columns and a plurality of rows, each of the plurality of 
columns representing the gray level value and each of the plurality of rows 
representing one of the plurality of slices in the scan, the entries corresponding to 
the number of times that the gray level occurs in the corresponding slice; 

determine a sum of tissue of voxels in the segmented lung image from the 
matrix entries; and 

multiply the sum of tissue of voxels in the segmented lung image by the 
volume of each voxel. 

20. An apparatus for measuring volume of tissue as defined by Claim 1 9, 
wherein the sum of tissue of voxels is calculated by summing the product of each 
matrix entry multiplied by the corresponding gray level value divided by a gray 
level value assigned for tissue. 

53 



wo 03/024184 



PCT/US02/29366 



21 . An article of manufacture for generating a lung mask for segmenting a lung 
image from voxel data containing the lung image, noise, solid components and 
surrounding background information, the article comprising: 

a machine readable medium containing one or more programs which when 
executed implement the steps of: 

(a) applying a median filter and a mean filter to the voxel data to reduce 

noise; 

(b) thresholding the noise reduced voxel data to identify solid 
components and other structures; 

(c) identifying the surrounding background in the noise reduced voxel 

data; 

(d) deleting the surrounding background from the noise reduced voxel 

data; 

(e) labeling coimected components of the voxel data; 

(f) determining the largest connected components to select a lung 
region from the voxel data, the lung region having a geometric form; and 

(g) morphological filtering the voxel data to refine the geometric form 
of the lung region. 

22. An article of manufacture for generating a lung mask as defined by Claim 
21, wherein the voxel data is associated with a plurality of slices through a 
patient's lung, the slices beginning at about the patient's shoulders; and 

step (a) is performed for the first 25 percent of the plurality of slices only 
whereby the computation time is substantially reduced. 

23. An article of manufacture for generating a lung mask as defined by Claim 
21 5 wherein the median filter has a size of 4x4. 

24. An article of manufacture for generating a lung mask as defined by Claim 
21, wherein the mean filter has a size of 1x3. 
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25. An article of manufacture for generating a lung mask as defined by Claim 
21, wherein the voxel data in step (b) is thresholded at a gray level of about 500. 

26. An article of manufacture for generating a lung mask defined by Claim 21, 
wherein the largest coimected components are associated with more than about 1 
percent of the voxel data. 

27. An article of manufacture for generating a lung mask as defined by Claim 
21, wherein the voxel data is associated with a plurality of slices, the slices being 
divided into a first end region, a middle region, and a second end region; and 

the morphological filtering being performed with 2D circular filter having a 
first diameter in the first end region and the second end region; and the 
morphological filtering being perfomied with 2D circular filter having a second 
diameter which is about twice the first diameter in the middle region. 

28. An article of manufacture for measuring lung volume from a segmented 
lung image obtained fi:om a scan which includes a plurality of slices of voxel data, 
the voxel data having a gray level value and a volume associated therewith, the 
article comprising: 

a machine readable medium containing one or more programs which when 
executed implement the steps of: 

generating a matrix of entries from the segmented lung image, the matrix 
having a plurality of columns and a plurality of rows, each of the plurality of 
columns representing the gray level value and each of the plurality of rows 
representing one of the plurality of slices in the scan, the entries corresponding to 
the number of times that the gray level occurs in the corresponding slice; 

determining a number of voxels in the segmented lung image from the 
matrix entries; and 

multiplying the number of voxels in the segmented lung image by the 
volume of each voxel. 
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29. An article of manufacture for measuring volume of tissue in a segmented 
Ixmg image obtained firom a scan including a plurality of slices of voxel data, the 
voxel data having a gray level value aad a volume, the article comprising: 

a machine readable medium containing one or more programs which when 
executed implement the steps of: 

generating a matrix of entries from the segmented lung image, the matrix 
having a plurality of columns and a plurality of rows, each of the plurality of 
columns representing the gray level value and each of the plurality of rows 
representing one of the plurality of slices in the scan, the entries corresponding to 
the number of times that the gray level occurs in the corresponding slice; 

determining a sum of tissue of voxels in the segmented lung image from the 
matrix entries; and 

multiplying the sum of tissue of voxels in the segmented Ixmg image by the 
volume of each voxel. 

30. An article of manufacture for measuring volume of tissue as defined by 
Claim 29, wherein tlie sum of tissue of voxels is calculated by summing the 
product of each matrix entry multiplied by the corresponding gray level value 
divided by a gray level value assigned for tissue. 

31. A method for finding the location, P\ and size, r\ of a pulmonary nodule in 
a high-resolution computed tomography (CT) image, the method comprising: 

(a) selecting a set of initial processing parameters including an initial 
location. Pi, an initial size, ri, and target value, T; 

(b) computing an initial new location/} of the nodule with a locator 

template ftmction; 

(c) incrementally increasing size, r/; 

(d) computing a new location of the nodule with the locator template 
ftmction; 

(e) computing a size metric with a sizing template ftmction; 
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(f) repeating steps (c) through (e) untU the size metric is less than the 
target value; 

(g) returning the location, P\ and the size, r\ from the previous 
iteration of steps (c) through (e). 

32. A method for finding the location, P, and size, r, of a pulmonary nodule in 
a high-resolution computed tomography image, the method comprising: 

(a) windowing the image to ignore bone structures; 

(b) selecting a locator template function and a sizing template function; 

(c) selecting a set of initial processing parameters including: an initial 
location, P\ size, r; and termination criteria; 

(d) performing a search to determine a maximum response of the 
locator template function; 

(e) determining a response of the sizing template function and 
comparing the response to the termination criteria; 

(f) incrementally increasing the size, r, only if the termination criteria 
has not been satisfied and repeating steps d and e; and 

(g) outputting the location, P, and size, r, of the nodule. 

33. A method for finding the location and size of a pulmonary nodule as 
defined in Claim 32, wherein the windowing comprises clipping the image at 
intensities over about 1000, 

34. A method for finding the location and size of a pulmonary nodule as 
defined in Claim 32, wherein the locator template function is a Gaussian template 
function. 

35. A method for finding the location and size of a pulmonary nodule as 
defined in Claim 32, wherein the locator template function is a Laplacian of the 
Gaussian template function. 
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36. A method for finding the location and size of a pulmonary nodule as 
defined in Claim 32, wherein the locator template function is a difference of 
Gaussians template function. 

37. A method for finding the location and size of a pulmonary nodule as 
defined in Claim 32, wherein the template function has at least four parameters 
corresponding to the x-location, y-location, z-location, and radius. 

38. A method for finding the location and size of a pulmonary nodule as 
defined in Claim 32, wherein the initial location, P, is calculated from the image. 

39. A method for finding the location and size of a pulmonary nodule as 
defined in Claim 32, wherein the initial location, is specified by a user. 

40. A method for finding Ihe location and size of a pulmonary nodule as 
defined in Claim 32, wherein the search is performed by a hill climbing method. 

41. A method for finding the location and size of a pulmonary nodule as 
defined in Claim 32, wherein the search is performed by Powell's method. 

42. A pulmonary nodule finding apparatus for finding the location, P\ and size, 
r\ of a pulmonary nodule in a high-resolution computed tomography (CT) image, 
the pulmonary nodule finding apparatus comprising: 

a nodule finding unit configured to: 

(a) select a set of initial processing parameters including an initial 
location. Pi, an initial size, ri, and target value, T; 

(b) compute an initial new location of the nodule with a locator 

template function; 

(c) incrementally increase size, r,-; 

(d) compute a new location P^ of the nodule with the locator template 
function; 

(e) compute a size metric with a sizing template function; 
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(f) repeat steps (c) through (e) until the size metric is less than the 
target value; 

(g) return the location, P\ and the size, r', from the previous iteration of 
steps (c) through (e). 

43. A pulmonary nodule jSnding apparatus for finding the location, P, and size, 
r, of a pulmonary nodule in a high-resolution computed tomography image, the 
pulmonary nodule finding apparatus comprising: 

a nodule finding unit configured to: 

(a) window the image to ignore bone structures; 

(b) select a locator template function and a sizing template function; 

(c) select a set of initial processing parameters including: an initial 
location, P; size, r; and termination criteria; 

(d) perform a search to determine a maximum response of the locator 
template function; 

(e) determine a response of the sizing template function and comparing 
the response to the termination criteria; 

(f) incrementally increase the size, r, only if the termination criteria has 
not been satisfied and repeating steps d and e; and 

(g) output the location, P, and size, r, of the nodule. 

44. A pulmonary nodule finding apparatus as defined in Claim 43, wherein the 
image is clipped at intensities over about 1000 to window the image. 

45. A pulmonary nodule finding apparatus as defined in Claim 43, wherein the 
locator template function is a Gaussian template function. 

46. A pulmonary nodule finding apparatus as defined in Claim 43, wherein the 
locator template function is a Laplacian of the Gaussian template function. 

47. A pulmonary nodule finding apparatus as defined in Claim 43, wherein the 
locator template function is a difference of Gaussians template function. 
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48. A pulmonary nodule finding apparatus as defined in Claim 43, wherein the 
template function has at least four parameters corresponding to the x-location, y- 
location, z-location, and radius, 

49. A pulmonary nodule finding apparatus as defined in Claim 43, wherein the 
initial location, P, is calculated from the image. 

50. A pulmonary nodule finding apparatus as defined in Claim 43 , wherein the 
initial location, P, is specified by a user. 

51. A pulmonary nodule finding apparatus as defined in Claim 43 , wherein the 
search is performed by a hill climbing method. 

52. A pulmonary nodule finding apparatus as defined in Claim 43, wherein the 
search is performed by Powell's method. 

53. An article of manufacture for finding the location, P\ and size, r\ of a 
pulmonary nodule in a high-resolution computed tomography (CT) image, the 
article comprising: 

a machine readable medium containing one or more programs which when 
executed implement the steps of: 

(a) selecting a set of initial processing parameters including an initial 
location, Pi, an initial size, ri, and target value, T; 

(b) computing an initial new location jFJ. of the nodule with a locator 

template function; 

(c) incrementally increasing size, r/; 

(d) computing a new location P, of the nodule with the locator template 
function; 

(e) computing a size metric with a sizing template function; 

(f) repeating steps (c) through (e) xmtil the size metric is less than the 
target value; 
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(g) returning the location, P\ and the size, r', from the previous 
iteration of steps (c) through (e). 

54. An article of manufacture for finding the location, P, and size, r, of a 
pulmonary nodule in a high-resolution computed tomography image, the article 
comprising: 

a machine readable medixmi containing one or more programs which when 
executed implement the steps of: 

(a) windowing the image to ignore bone structures; 

(b) selecting a locator template function and a sizing template function; 

(c) selecting a set of initial processing parameters including: an initial 
location, P; size, r; and termination criteria; 

(d) performing a search to determine a maximum response of the 
locator template function; 

(e) determining a response of the sizing template function and 
comparing the response to the termination criteria; 

(f) incrementally increasing the size, r, only if the termination criteria 
has not been satisfied and repeating steps d and e; and 

(g) outputting the location, P, and size, r, of the nodule. 

55. An article of manufacture for finding the location and size of a pulmonary 
nodule as defined in Claim 54, wherein the windowing comprises clipping the 
image at intensities over about 1000. 

56. An article of manufacture for finding the location and size of a pulmonary 
nodule as defined in Claim 54, wherein the locator template function is a Gaussian 
template function. 

57. An article of manufacture for finding the location and size of a pulmonary 
nodule as defined in Claim 54, wherein the locator template function is a Laplacian 
of the Gaussian template function. 
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58. An article of manufacture for finding the location and size of a pulmonary 
nodule as defined in Claim 54, wherein the locator template function is a 
difference of Gaussians template function. 

59. An article of manufacture for finding the location and size of a pulmonary 
nodule as defined in Claim 54, wherein the template function has at least four 
parameters corresponding to the x-location, y-location, z-location, and radius. 

60. An article of manufacture for finding the location and size of a pulmonary 
nodule as defined in Claim 54, wherein the initial location, is calculated fi:om 
the image. 

61. An article of manufacture for finding the location and size of a pulmonary 
nodule as defined in Claim 54, wherein the initial location, P, is specified by a 
user. 

62. An article of manufacture for finding the location and size of a pulmonary 
nodule as defined in Claim 54, wherein the search is performed by a hill climbing 
method. 

63. An article of manufacture finding the location and size of a pulmonary 
nodule as defined in Claim 54, wherein the search is performed by PowelFs 
method. 

64. A method for registering 3-d images of a pulmonary nodule from a high- 
resolution computed tomography (CT) scans, the images being in a floating point 
pixel-format associated with a 6-dimensional parameter space and including a first 
image {im\) obtained at time-1 and a second image (/mi) obtained at time-2, the 
method comprising the steps of: 

(a) calculating initial rigid-body transformation parameters for a rigid- 
body transformation on the first image (/mi); 
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(b) determining the optimum rigid-body transformation parameters by 
calculating a registration metric between the second image (imi) and the rigid-body 
transformation on the first image (im\)\ and 

(c) generating a registered image from the optimum rigid-body 
transformation parameters. 

65. A method as defined in Claim 64, wherein step (a) is preceded by masking 
one of the images by setting pixels to a background value. 

66. A method as defined in Claim 65, wherein the background value is about - 
1000. 

67. A method as defined in Claim 64, wherein the registration metric is 
minimized. 

68. A method as defined in Claim 64, wherein the registration metric is 
maximized. 

69. A method as defined in Claim 64, wherein the registration metric is 
calculated by 

transforming the first image {im{) with the initial rigid-body 
transformation parameters to obtain a transformed first image (?>?ziO; 

calculating the registration metric as a correlation (C) between the 
transformed first image {im\^) and the second image (z>?Z2); and 

searching for the maximum correlation (C) in the 6-dimensional 
parameter space. 

70. A method as defined in Claim 64, wherein the registration metric is 
calculated by 

transforming the first image {im\) with the initial rigid-body 
transformation parameters to obtain a transformed first image {imuy, 

calculating the registration metric as a mean-squared-difference 
(MSD) between the transformed first image {im\i) and the second image (zm2); and 
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searching for the minimmn mean-squared-difference (MSD) in the 
6-diniensional parameter space. 



71 . A method as defined in Claim 69, wherein the transforming of the first 
image (imi) to obtain the transformed first image (imu) is a mapping of a point v 
in 3-d space to a point v* in transformed space defined by : 



wherein Rj^, Ry^ and R^ are rotation matrices defined as: 



1 0 0 
0 cos(r^) -sin(r^) 
0 sin(rj cos(rJ 



cos(r^) 0 sin(r^) 

0 1 0 
-sin(r^) 0 cosCr^) 



cos(r^) 
sin(rj 
0 



-sin(rj 0 
cos(r^) 0 
0 1 



72. A method as defined in Claim 70, wherein the transforming of the first 
image (zmi) to obtain the transformed first image (imu) is a mapping of a point v 
in 3-d space to a point v* in transformed space defined by : 



wherein R-^, Ry, and R^ are rotation matrices defined as: 
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10 0 
0 cos(r^) -sm(r^) 
0 sin(r^) cos(r^) 



cos(r^) 0 sm(r^) 

0 1 0 
-sin(ry) 0 cos(ry) 



cos(rJ -sin(rj 0 
i?, = sin(rj cos(;;) 0 
0 0 1 



73. A method as defined in Claim 64, wherein initial rigid-body transformation 
parameters include six parameters (tx,tyjz,rx^fy,rz) respectively defined as 
translation in x, translation in y, translation in z, rotation about the x-axis, rotation 
about the y-axis, and rotation about the z-axis; 

wherein the initial rotation parameters {rx,}y,7'z) are all set to zero; and the 
initial translation parameters (tx,ty,tz,) are set so that the nodule in the first image 
(imi) overlaps the nodule in the second image (im2) during the initial calculation of 
the registration metric. 

74. A method as defined in Claim 73, wherein the initial translation parameters 
(tx,ty,tz,) are set to a difference between the center of the first image (im\) and the 
center of the second image (imz). 

75. A method as defined in Claim 73, wherein the initial translation parameters 
(pc^ty^tz) are set to a difference between the center of mass of the first image (imi) 
and the center of mass of the second image (im2)* 



65 



wo 03/024184 



PCT/US02/29366 



76. A method as defined in Claim 69, wherein the searching is conducted by 
calculating the correlation (C) for every possible set of rigid-body transformation 
parameters. 

77. A method as defined in Claim 70, wherein the searching is conducted by 
calculating the mean-squared-difference (MSD) for every possible set of rigid- 
body transformation parameters. 

78. A method as defined in Claim 69, wherein the searching is conducted by a 
Hill-Climbing search method. 

79. A method as defined in Claim 70, wherein the searching is conducted by a 
Hill-Climbing search method. 

80. A method as defined in Claim 69, wherein the searching is conducted by 
PowelPs method. 

81 . A method as defined in Claim 70, wherein the searching is conducted by 
Powell's method. 

82. A registering apparatus for registering 3-d images of a pulmonary nodule 
from a high-resolution computed tomography (CT) scans, the images being in a 
floating point pixel-format associated with a 6-dimensional parameter space and 
including a first image (im{) obtained at time-1 and a second image (2/^2) obtained 
at time-2, the registering apparatus comprising: 

a registering unit configured to: 

(a) calculate initial rigid-body transformation parameters for a rigid- 
body transformation on the first image (im\); 

(b) determine the optimum rigid-body transformation parameters by 
calculating a registration metric between the second image (im2) and the rigid-body 
transformation on the first image (im\); and 
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(c) generate a registered image from the optimimi rigid-body 
transformation parameters. 

83. A registering apparatus as defined in Claim 82, wherein one of the images 
is initially masked by setting pixels to a background value. 

84. A registering apparatus as defined in Claim 83, wherein the background 
value is about -1000. 

85. A registering apparatus as defined in Claim 82, wherein the registration 
metric is minimized. 

86. A registering apparatus as defined in Claim 82, wherein the registration 
metric is maximized. 

87. A registering apparatus as defined in Claim 82, wherein the registration 
metric is calculated by 

transforming the first image (imi) with the initial rigid-body 
transformation parameters to obtain a transformed first image (imu); 

calculating the registration metric as a correlation (C) between the 
transformed first image (iniu) and the second image (im2); and 

searching for the maximum correlation (C) in the 6-dimensional 
parameter space. 

88. A registering apparatus as defined in Claim 82, wherein the registration 
metric is calculated by 

transforming the first image (imi) with the initial rigid-body 
transformation parameters to obtain a transformed first image (mtu); 

calculating the registration metric as a mean-squared-difference 
(MSD) between the transformed first image (z>72u) and the second image (2^2); and 

searching for the minimum mean-squared-difference (MSD) in the 
6-dimensional parameter space. 
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89. A registering apparatus as defined in Claim 87, wherein the transforming of 
the first image (imi) to obtain the transformed first image (iniu) is a mapping of a 
point V in 3-d space to a point v' in transformed space defined by : 

wherein R^, Ry, and i?^ are rotation matrices defined as: 




10 0 ' 

0 cos(r^) -sin(r^) 
0 sin(r,) cos(r^) 



cos(ry) 0 sin(ry) 

0 1 0 
-sin(r^) 0 cos(r^) 



R, = 



"cos(7;) 
sin(rj 
0 



— s\n{r^) 0 
cos(rJ 0 
0 1 



90. A registering apparatus as defined in Claim 88, wherein the transforming of 
the first image (zmi) to obtain the transformed first image {imu) is a mapping of a 
point V in 3-d space to a point v' in transformed space defined by : 



'yRz 



V+ 



wherein R^^ Ry, and i?^ are rotation matrices defined as: 



10 0 

0 cos(rJ -sin(r,) 
0 sin(r,) cos(r^) 
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cos(r^) 0 sm(r^,) 

0 1 0 
~sin(r^) 0 cos(r^) 



cos(r^ ) 
siii(rj 
0 



cos(r^) 0 
0 1 



91. A registering apparatus as defined in Claim 82, wherein initial rigid-body 
transformation parameters include six parameters (tXyty^tz^rx,ry,rz) respectively 
defined as translation in x, translation in y, translation in z, rotation about the x- 
axis, rotation about the y-axis, and rotation about the z-axis; 

wherein the initial rotation parameters {rxy7y,rz) are all set to zero; and the 
initial translation parameters (tx^ty.tz^) are set so that the nodule in the first image 
(imi) overlaps the nodule in the second image (imz) during the initial calculation of 
the registration metric, 

92. A registering apparatus as defined in Claim 91, vv^herein the initial 
translation parameters (tx,ty,tz,) are set to a difference between the center of the 
first image (/mi) and the center of the second image (rm2). 

93. A registering apparatus as defined in Claim 91 , wherein the initial 
translation parameters (tx^tyytz) are set to a difference between the center of mass of 
the first image (imi) and the center of mass of tiie second image (im2), 

94. A registering apparatus as defined in Claim 87, wherein the searching is 
conducted by calculating the correlation (C) for every possible set of rigid-body 
transformation parameters. 

95. A registering apparatus as defined in Claim 88, wherein the searching is 
conducted by calculating the mean-squared-difference (MSD) for every possible 
set of rigid-body transformation parameters. 

69 



wo 03/024184 



PCT/US02/29366 



96. A registering apparatus as defined in Claim 87, wherein the searching is 
conducted by a Hill-Climbing search method. 

97. A registering apparatus as defined in Claim 88, wherein the searching is 
conducted by a Hill-Climbing search method. 

98. A registering apparatus as defined in Claim 87, wherein the searching is 
conducted by Powell's method. 

99* A registering apparatus as defined in Claim 88, wherein the searching is 
conducted by Powell's method. 

100. An article of manufacture for registering 3-d images of a pulmonary nodule 
firom a high-resolution computed tomography (CT) scans, the images being in a 
floating point pixel-format associated with a 6-dimensional parameter space and 
including a first image (imi) obtained at time-1 and a second image (iniz) obtained 
at time-2, the article comprising: 

a machine readable medium containing one or more programs which when 
executed implement the steps of: 

(a) calculating initial rigid-body transformation parameters for a rigid- 
body transformation on the first image (imi); 

(b) detemiining the optimum rigid-body transformation parameters by 
calculating a registration metric between the second image (miz) and the rigid-body 
transformation on the first image (imi); and 

(c) generating a registered image firom the optimum rigid-body 
transformation parameters. 

101 - An article of manufacture for registering 3-d images of a pulmonary nodule 
as defined in Claim 100, wherein step (a) is preceded by masking one of the 
images by setting pixels to a background value. 
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102. An article of manufacture for registering 3-d images of a pulmonary nodule 
as defined in Claim 101, wherein the background value is about -1000. 

103. An article of manufacture for registering 3-d images of a pulmonary nodule 
as defined in Claim 100, wherein the registration metric is minimized. 

104. An article of manufacture for registering 3-d images of a pulmonary nodule 
as defined in Claim 100, wherein the registration metric is maximized. 

105. An article of manufacture for registering 3-d images of a pulmonary nodule 
as defined in Claim 100, wherein the registration metric is calculated by 

transforming the first image {im\) with the initial rigid-body 
transformation parameters to obtain a transformed first image (zmu); 

calculating the registration metric as a correlation (C) between the 
transformed first image {im\i) and the second image {imi); and 

searching for the maximum correlation (C) in the 6-dimensional 
parameter space. 

106. An article of manufacture for registering 3-d images of a pulmonary nodule 
as defined in Claim 100, wherein the registration metric is calculated by 

transforming the first image (/mi) with the initial rigid-body 
transformation parameters to obtain a transformed first image {imu)\ 

calculating the registration metric as a mean-squared-difference 
(MSD) between the transformed first image (iniu) and the second image (iniT); and 

searching for the minimum mean-squared-difference (MSD) in the 
6-dimensional parameter space. 

107. An article of manufacture for registering 3-d images of a pulmonary nodule 
as defined in Claim 105, wherein the transforming of the first image (zmi) to obtain 
the transformed first image {im\^ is a mapping of a point v in 3-d space to a point 
V in transformed space defined by : 
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wherein R^^ Ry, and R^ are rotation matrices defined as: 



1 0 0 

0 cos(r^) -sin(r^) 
0 sin(r^) cos(r^) 



R. 



cos(r^) 0 sin(r^) 

0 1 0 
-sin(r^) 0 cos(r^) 



R. = 



cos(;;) -sin(/;) 0 
sin(r^) cos(/;) 0 
0 0 1 



108. An article of manufacture for registering 3-d images of a pulmonary nodule 
as defined in Claim 106, wherein the transforming of the first image (mi\) to obtain 
the transformed first image (imu) is a mapping of a point v in 3-d space to a point 
v' in transformed space defined by : 



V=RJtyRzV+ 



wherein i?^, Ry, and R^ are rotation matrices defined as: 



"l 0 0 
if^ = 0 cos(r^) -sin(r^) 
0 sin(r^) cos(r^) 
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cos(r^) 0 sin(r^) 

0 1 0 
-sin(r^) 0 cos(r^)_ 



"cos(rJ 
sin(rj 
0 



-sin(rj 0 
cos(rJ 0 
0 1 



109. An article of manufacture for registering 3-d images of a pulmonary nodule 
as defined in Claim 100, wherein initial rigid-body transfomiation parameters 
include six parameters (tx,ty,tz,rx^ry^T''z) respectively defined as translation in x, 
translation in y, translation in z, rotation about the x-axis, rotation about the y-axis, 
and rotation about the z-axis; 

wherein the initial rotation parameters (rx,ry,rz) are all set to zero; and the 
initial translation parameters (tx,ty,tz,) are set so that the nodule in the first image 
(imi) overlaps the nodule in the second image (im2) during the initial calculation of 
the registration metric. 

110. An article of manufacture for registering 3-d images of a pulmonary nodule 
as defined in Claim 109, wherein the initial translation parameters (tx,ty,tz,) are set 
to a difference between the center of the first image (zmi) and the center of the 
second image (im2). 

111. An article of manufacture for registering 3-d images of a pulmonary nodule 
as defined in Claim 109, wherein the initial translation parameters (tx.tyjz) are set 
to a difference between the center of mass of the first image (imi) and the center of 
mass of the second image (2/^2). 

112. An article of manufacture for registering 3-d images of a pulmonary nodule 
as defined in Claim 105, wherein the searching is conducted by calculating the 
correlation (C) for every possible set of rigid-body transformation parameters. 
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113. An article of manufacture for registering 3-d images of a pulmonary nodule 
as dejBned in Claim 106, wherein the searching is conducted by calculating the 
mean-squared-difference (MSD) for every possible set of rigid-body 
transformation parameters. 

1 14. An article of manufacture for registering 3-d images of a pulmonary nodule 
as defined in Claim 105, wherein the searching is conducted by a Hill-Climbing 
search method. 

115. An article of manufacture for registering 3-d images of a pulmonary nodule 
as defined in Claim 106, wherein the searching is conducted by a Hill-Climbing 
search method, 

116. An article of maaufacture for registering 3-d images of a pulmonary nodule 
as defined in Claim 105, wherein the searching is conducted by Powell's method. 

117. An article of manufacture for registering 3-d images of a pulmonary nodule 
as defined in Claim 106, wherein the searching is conducted by Powell's method. 

118. A method for removing extraneous matter from an image, the image 
including a juxtapleural nodule, the method comprising the steps of: 

(a) providing an initial location P'; 

(b) calculating a spherical volume centered at the initial location P\ the 
spherical volume fitting inside the image; 

(c) calculating a center of mass COM of the spherical volume; 

(d) determining an initial direction d\ the initial direction rf' being 
directed towards the extraneous matter in accordance with the following equation 



^. COM-F 



CQM-P 



(e) initializing a current location P,- to be equal to the initial location P\ 
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(f) initializing a current direction di to be equal to the initial direction 

d\ 

(g) initializing a maximum ratio 7max? step size prior mass massi.u 
and prior change in mass A,_i; 

(h) moving the current location Pi by the step size s in the current 
direction di\ 

(i) determining an equation defining a plane A, the plane A being 
normal to the current direction du the plane A passing through the current location 
Pi\ 

(j) calculating a current mass massi of the nodule on a side of the plane 
A opposing that of the extraneous matter; 

(k) calculating a current change in mass A/ by subtracting the prior 
mass massi-i from the current mass massi\ 

(1) calculating a current ratio 7 in accordance with the following 
equation 



(m) setting the prior mass massi^\ equal to the current mass massi\ 
(n) setting the prior change in mass A/.i equal to the current change in 
mass A/; 

(o) comparing the current ratio 7 to the maximum ratio 7niax; 

(p) modifying the current direction di to minimize the current mass 
massu and performing steps (h)-(o) while the current ratio 7 is one of less than and 
equal to the maximum ratio 7max; and 

(q) modifying the current direction di^ performing steps (h)-(o)5 and 
outputting the area of the nodule partitioned by the plane A in response to the 
current ratio 7 being greater than the maximum ratio 7max. 
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119. A method for removing extraneous matter from an image, the image 
including a juxtapleural nodule as dejBned in Claim 118, wherein the extraneous 
matter includes a pleural surface. 

120. A method for removing extraneous matter from an image, the image 
including a juxtapleural nodule as defined in Claim 118, wherein step (g) fiirther 
includes the steps of: 

initializing the maximum ratio 7max to 0.5; and 
initializing the step size 5 to 1.5. 

121 . A method for removing extraneous matter from an image, the image 
including a juxtapleural nodule as defined in Claim 118, wherein step (j) fiirther 
includes the steps of: 

(r) defining the current location Pi as being visited; 
(s) determining on which side of the plane A the current location P/ is 
located; 

(t) terminating in reponse to the current location Pi not being located 
on a side of the plane opposing that of the extraneous matter; 

(u) defining the current location Pi as being part of a region of interest 
in response to the current location Pi being located on the side of the plane 
opposing that of the extraneous matter; and 

(v) performing steps (r)-(u) recursively using a location corresponding 
to at least one of six (6) one-pixel moves from the current location Pi, 

122. A method for removing extraneous matter from an image, the image 
including a juxtapleural nodule as defined in Claim 118, wherein step (p) fiirther 
includes the steps of: 

calculating recursively the current mass niassi of the nodule on a side of the 
plane A opposing that of the extraneous matter using at least one of six (6) 
directions and a step size si from the current location Pc, and 

defining the current direction di equal to the direction yielding the largest 
decrease in the current mass massf. 
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123 . A method for removing extraneous matter from an image, the image 
including a juxtapleural nodule as defined in Claim 118, wherein step (q) further 
includes the steps of: 

calculating recursively the current mass massi of the nodule on a side of the 
plane A opposing that of the extraneous matter using at least one of six (6) 
directions and a step size si from the current location Pi\ and 

defining the current direction di equal to the direction yielding the largest 
decrease in the current mass massi, 

124. A method for removing extraneous matter from an image, the image 
including a juxtapleural nodule as defined in Claim 118, wherein the initial 
location P' is located near a center of the nodule. 

125. A recursive apparatus for removing extraneous matter from an image, the 
image including a juxtapleural nodule, the recursive apparatus comprising: 

a processing unit, the processing unit being configured to: 

(a) accept an initial location P 

(b) calculate a spherical volume centered at the initial location P\ the 
spherical volume fitting inside the image; 

(c) calculate a center of mass COM of the spherical volume; 

(d) determine an initial direction d\ the initial direction d* being 
directed towards the extraneous matter in accordance with the following equation 

COM-F 



COM-P 



(e) initialize a current location Pi to be equal to the initial location P'\ 

(f) initialize a current direction dt to be equal to the initial direction J'; 

(g) initialize a maximum ratio 7maxj step size s, prior mass massi-u and 
prior change in mass Am; 
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(h) move the current location P/ by the step size s in the current 
direction df, 

(i) determine an equation defining a plane A, the plane A being normal 
to the current direction di, the plane A passing through the current location P,-; 

(j) calculate a current mass massi of the nodule on a side of the plane A 
opposing that of the extraneous matter; 

(k) calculate a current change in mass A/ by subtracting the prior mass 
jnassi-i from the current mass massi\ 

(1) calculate a current ratio 7 in accordance with tfie following equation 



(m) set the prior mass massi.\ equal to the current mass massi\ 

(n) set the prior change in mass Am equal to the current change in mass 



A,; 



(o) compare the current ratio 7 to the maximum ratio 7inax; 

(p) modify the current direction di to minimize the current mass massi^ 
and performing steps (h)-(o) while the current ratio 7 is one of less than and equal 
to the maximum ratio 7max; and 

(q) modify the current direction J/, performing steps (h)-(o), and 
outputting the area of the nodule partitioned by the plane A in response to the 
current ratio 7 being greater than the maximum ratio 7max. 

126. A recursive apparatus for removing extraneous matter from an image as 
defined in Claim 125, wherein the extraneous matter includes a pleural surface. 

127. A recursive apparatus for removing extraneous matter from an image as 
defined in Claim 125, wherein in step (g) the processing unit is configured to: 

initialize the maximum ratio 7max to 0.5; and 
initialize the step size to 1.5. 
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128. A recursive apparatus for removing extraneous matter from an image as 
defined in Claim 125, wherein in step (j) the processing imit is configured to: 

(r) define the current location Pi as being visited; 
(s) determine on which side of the plane A the current location Pi is 
located; 

(t) terminate in reponse to the current location P,- not being located on a 
side of the plane opposing that of the extraneous matter; 

(u) define the current location Pi as being part of a region of interest in 
response to the current location Pi being located on the side of the plane opposing 
that of the extraneous matter; and 

(v) perform steps (r)-(u) recursively using a location corresponding to 
at least one of six (6) one-pixel moves from the current location P/. 

129. A recursive apparatus for removing extraneous matter from an image as 
defined in Claim 125, wherein in step (p) the processing unit is configured to: 

calculate recursively the current mass massi of the nodule on a side of the 
plane A opposing that of the extraneous matter using at least one of six (6) 
directions and a step size s from the current location P,- 1; and 

define the current direction di equal to the direction yielding the largest 
decrease in the current mass masSi, 

130. A recursive apparatus for removing extraneous matter from an image as 
defined in Claim 125, wherein in step (q) the processing unit is configured to: 

calculate recursively the current mass massi of the nodule on a side of the 
plane A opposing that of the extraneous matter using at least one of six (6) 
directions and a step size Si from the current location Pi\ and 

defining the current direction di equal to the direction yielding the largest 
decrease in the current mass masSi, 

131. A recursive apparatus for removing extraneous matter from an image as 
defined in Claim 125, wherein the initial location P' is located near a center of the 
nodule. 
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132. An article of manufacture for removing extraneous matter from an image, 
the image including a juxtapleural nodule, the article comprising: 

a machine readable medium containing at least one program which when 
executed implements the steps of: 

(a) accepting an initial location P 

(b) calculating a spherical volvime centered at the initial location P', the 
spherical volume fitting inside the image; 

(c) calculating a center of mass COM of the spherical volume; 

(d) determining an initial direction d\ the initial direction rf' being 
directed towards the extraneous matter in accordance with the following equation 

COM-P' ^ 
pOM-Pf 

(e) initializing a current location Pi to be equal to the initial location P\ 

(f) initializing a current direction di to be equal to the initial direction 

d\ 

(g) initializing a maximum ratio Ymax^ step size prior mass inassuu 
and prior change in mass A/.i; 

(h) moving the current location Pf by the step size s in the current 
direction di\ 

(i) determining an equation defining a plane A, the plane A being 
normal to the current direction d^ the plane A passing through the current location 
Pi\ 

(j) calculating a current mass massi of the nodule on a side of the plane 
A opposing that of the extraneous matter; 

(k) calculating a current change in mass A/ by subtracting the prior 
mass massi.\ from the current mass rnassi; 

(1) calculating a current ratio 7 in accordance with the following 
equation 
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(m) setting the prior mass massi.i equal to the current mass massi\ 
(n) setting the prior change in mass A/_i equal to the current change in 
mass A/; 

(o) comparing the current ratio 7 to the maximum ratio 7max; 

(p) modifying the current direction dt to minimize the current mass 
massu and performing steps (h)-(o) while the current ratio 7 is one of less than and 
equal to the maximum ratio 7max; and 

(q) modifying the current direction di, performing steps (h)-(o), and 
outputting the area of the nodule partitioned by the plane A in response to the 
current ratio 7 being greater than the maximum ratio 7max. 

133. An article of manufacture for removing extraneous matter from an image as 
defined in Claim 132, wherein the extraneous matter includes a pleural surface. 

134. An article of manufacture for removing extraneous matter firom an image as 
defined in Claim 132, wherein in step (g) the article further implements the steps 
of: 

initializing the maximum ratio 7niax to 0,5; and 
initializing the step size s to 1.5. 

135. An article of manufacture for removing extraneous matter firom an image as 
defined in Claim 132, wherein in step (j) the article further implements the steps 
of: 

(r) defining the current location P/ as being visited; 
(s) determining on which side of the plane A the current location Pi is 
located; 

(t) terminating in reponse to the current location Pi not being located 
on a side of the plane opposing that of the extraneous matter; 
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(u) defining the current location Pi as being part of a region of interest 
in response to the current location Pi being located on the side of the plane 
opposing that of the extraneous matter; and 

(v) performing steps (r)-(u) recursively using a location corresponding 
to at least one of six (6) one-pixel moves from the current location P/. 

136. An article of manufacture for removing extraneous matter from an image as 
defined in Claim 132, wherein in step (p) the article ftirther implements the steps 
of: 

calculating recursively the current mass massi of the nodule on a side of the 
plane A opposing that of the extraneous matter using at least one of six (6) 
directions and a step size si from the current location P/; and 

defining the current direction di equal to the direction yielding the largest 
decrease in the current mass masSi. 

137. An article of manufacture for removing extraneous matter from an image as 
defined in Claim 132, wherein in step (q) the article fiirther implements the steps 
of: 

calculating recursively the current mass masSi of the nodule on a side of the 
plane A opposing that of the extraneous matter using at least one of six (6) 
directions and a step size Si from the current location P/; and 

defining the current direction di equal to the direction yielding the largest 
decrease in the current mass massi, 

138. An article of manufacture for removing extraneous matter from an image as 
defined in Claim 132, wherein the initial location P' is located near a center of the 
nodule. 
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FIG. 5 



Lung Mask Generation Algorithm 

1. Apply median and mean filters to the original scan 

to reduce noise. 
2. Segment the scan into solid components and other 

structures via thresholding. 
3. Identify and delete the surrounding background. 

(i.e. region that is not part of the thorax). 

4. Label all connected components and select the lung 
region (the largest component). 

5. Refine the geometric form of the resulting region 
by morphological filtering. 
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FIG. 6 
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FIG. 15 



Algorithm for finding the location and size of a nodule 

l.Initial pointPl = (xpjj, z^) and size r^. 

2. Increase radius: = r ^-^j + 

3. Compute best location: P search (P^. j 

4. Compute match (size metric): jj = size (P^ r) 

5. Repeat 2 tlirough 4 until 7/ < T 

6. P' = P .,r'-r^. 
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Detailed Algorithm for finding the location P and size r of a nodule 

1. Given a high-resolution 3d CT image of a nodule. 

2. Window the image to ignore any bone structures. 

3. Define the locator and sizing template functions. 

4. Define the initial location, P, 

5. Define a termination criteria: target^ ^djud a, 

6. Set rstaj't = 1, rstep = 1, rniax, 7 = oo. 

7. Wliile ( I ^-target \ > e) and {rstep > a) 

(a) For rf=rstart to rmax with increment of rstep 

i. Set the locator template radius to rj. 

ii. Calculate the location: Pi = search(Pi-\^rf) 

iii. Set the sizing template radius to rf and the location to Pj. 

iv. Calculate size: 7 — the response of the sizing template function. 
V. If (7 < target)^ then exit the for loop (go to 7b). 

(b) Set rstep = rstep I a 

(c) Set rstart = r/-l + rstep 

(d) SetP/ -=Pz-l 

8. Output the location and radius of the nodule, Pf, rf . 
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FIG. 23 




wo 03/024184 



PCT/US02/29366 



19/26 



FIG. 24 
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Hill-Climbing Search Algorithm for the nodule finding Algorithm 

1. Given an initial condition, initial and a stepsize, dt, 

2. par am = initial 

3. Evaluate the locator template function for movement of dt 
and - dt in the three directions from param, 

4. If a smaller response is calculated, then move param in the 
appropriate direction for the largest decrease and goto 3. 

5. Return the minimizing parameters, /?araw. 
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Algorithm for registering a first image, im\, and a second image, 
z>W25 using a rigid-body transformation 

1. Convert both input images, im\ and imi, to floating point 
pixels. 

2. If desired, mask either image by setting pixels to a 
background value. 

3. Calculate the initial conditions for the rigid-body 
transformation on imi. 

4. Minimize (or maximize) the metric between im2 and the 
rigid-body transformation of imi by changing the 
transformation parameters. 

5. Using the best rigid-body transformation parameters, 
transform imi, convert it back to the original pixel format, and 
output the image. 
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Hill-Climbing Search Algorithm for determining the best rigid-body transformation 
parameters 

1. Given an initial condition, initial, the translation stepsize, dt, the rotation 

stepsize, dr, and the number of passes, p 

2. Inflate the stepsizes: dt = dt "^IP"^, dr ^^dr"^ IP"^ 

3. par am == initial 

4. For a == 1 to p 

(a) Evaluate the metric for movement from param in 
the three translation directions with steps dt and - 
dt pixels. 

(b) Evaluate the metric for movement from param in 
the three rotation directions with steps dr and - dr 
degrees. 

(c) If a smaller metric value is calculated in (a) or (b), 
then move param in the appropriate direction for 
the largest decrease and goto (a). 

(d) Update stepsizes: dt - dt / 2, dr - dr / 2 

5. Return the minimizing parameters, param. 
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Algorithm for removing the pleural-surface given a location P' near the center of the 
nodule. 

1 . Calculate the center of mass, COM, of the largest spherical volume centered at P' 
that fits inside the image, COM P' 

2. Determine the initial direction towards the pleural wall: <^ — ' ^^qq^ — PTT 

3. Initialize Pq -=^P'\d^d' 

4. Set 7 j^^^j^ = 0.5; 7 = 0; step = L5; 

5. Wliile'y<7^;^ 

(a) Move location: Pf=Pi_Y\-step^d 

(b) Calculate the plane A normal to d and passing through point Pi 

(c) mass I = size of the cut nodule (the connected component region containiag 
point Pi and located on the negative side of plane Ai) 

(d) Calculate the change in mass: H^f^massi -massi^i 

A; 

(e) Calculate the ratio change: 7 = -r— ^ — 1 

(f) if 7 > jj^^ , then change d to reorientate the plane Ai to minimize the size of 
the cut nodule. Recalculate massp Ai^ and 7 . 

6. Retum the cut nodule using plane Ai^i. 
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FIG. 30E 



FIG. 30F 
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FIG. 31 



Algorithm for Recursively finding the cut nodule region, given the starting point 
p — P\ and the plane and the image im 

1 . Mark the pixel p as visited. 

2. If point p is not behind plane A then exit. 

3. Mark pixel located at p as part of the region. 

4. For each of the six possible one-pixel moves, calculate a new position p\ If the 
pixel p'is foreground and has not been visited, then recursively call the algorithm 
using/? =jt7'. 



Hill-Climbing Search Algorithm for removing the pleural surface 

1. Given the initial direction, and a stepsize, step^. 

2. cIq 

3. For each of tlie six possible stepd steps fi'om d, calculate the size of the cut-nodule. 

4. If a smaller size is calculated, then move d in the appropriate direction for the 

largest decrease, normalize d, and goto 3. 

5. Return the minimizing direction d. 




